Télécharger rigix.eso

Retour à la liste

Numérotation des lignes :

  1. C RIGIX SOURCE GG250959 17/12/21 21:15:16 9680
  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.  
  482. XGAU = 0.D0
  483. YGAU = 0.D0
  484. ZGAU = 0.D0
  485. C
  486. i6zz = 3
  487. IF (IDIM.EQ.3) i6zz = 4
  488. C
  489. c do ienr7=1,NBENRMAX
  490. do ienr7=1,NBENRJ
  491. do inod7=1,NBNN1
  492. c do i6=1,6
  493. CYT do i6=1,3
  494. do i6=1,i6zz
  495. LV7WRK(ienr7,1,i6,inod7) = 0.D0
  496. LV7WRK(ienr7,2,i6,inod7) = 0.D0
  497. enddo
  498. enddo
  499. enddo
  500.  
  501.  
  502. c*********************************************************
  503. c Calcul des fonction de forme std dans repere local
  504. c*********************************************************
  505.  
  506. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  507. c (et donc sur les Ni std)
  508. DO 2510 I=1,NBNN1
  509.  
  510. C++++ Calcul des Ni std
  511. c (rappel: 1:Ni 2:Ni,qsi 3:Ni,eta avec i=1,4)
  512. SHPWRK(1,I) = SHPTOT(1,I,KGAU)
  513. SHPWRK(2,I) = SHPTOT(2,I,KGAU)
  514. SHPWRK(3,I) = SHPTOT(3,I,KGAU)
  515. IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)
  516.  
  517. C++++ Calcul des coordonnees dans le repere global
  518. XGAU = XGAU + (SHPWRK(1,I)*XE(1,I))
  519. YGAU = YGAU + (SHPWRK(1,I)*XE(2,I))
  520. IF (IDIM.EQ.3) ZGAU = ZGAU + (SHPWRK(1,I)*XE(3,I))
  521.  
  522. 2510 CONTINUE
  523. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  524.  
  525.  
  526.  
  527. c*********************************************************
  528. c Passage des fonctions de forme std dans repere global
  529. c*********************************************************
  530.  
  531. C++++ CALCUL DES Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  532. CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  533. c if(J.eq.1.and.KGAU.eq.1)
  534. c &write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1)
  535.  
  536. c*********************************************************
  537. c Si on est pas du tout enrichi on peut sauter une grosse partie
  538. c*********************************************************
  539. if(NBENRJ.eq.0) goto 2999
  540.  
  541. c*********************************************************
  542. c Calcul des level set + leurs derivees dans repere global
  543. c*********************************************************
  544.  
  545. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  546. do 2520 ienr=1,NBENRJ
  547.  
  548. c MELVA1=MCHAM1.IELVAL(IENR)
  549. c segact,MELVA1
  550.  
  551. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  552. DO 2521 I=1,NBNN1
  553.  
  554. C++++ Le I eme noeud est-il ienr-enrichi?
  555. mlree1= TLREEL(IENR,I)
  556. if(mlree1.eq.0) goto 2521
  557.  
  558.  
  559. c Calcul du repere local de fissure(=PSI,PHI)
  560. c (rappel: 1,1:psi 1,2:phi 2,1 psi,x ...etc...)
  561. do 2522 inode=1,NBNN1
  562. c pour le H-enrichissement, on n a pas gardé PSI (inutile)
  563. if (ienr.ne.1) then
  564. c valeur de PSI au inode^ieme noeud
  565. xpsi1 = mlree1.PROG(inode)
  566. c qu on multiplie par la valeur de Ni^std au pt de G considéré
  567. LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
  568. & + (SHPWRK(1,inode)*xpsi1)
  569. LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
  570. & + (SHPWRK(2,inode)*xpsi1)
  571. LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
  572. & + (SHPWRK(3,inode)*xpsi1)
  573. IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  574. & + (SHPWRK(4,inode)*xpsi1)
  575. c valeur de PHI au inode^ieme noeud
  576. xphi1 = mlree1.PROG(NBNN1+inode)
  577. else
  578. xphi1 = mlree1.PROG(inode)
  579. endif
  580. LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
  581. & + (SHPWRK(1,inode)*xphi1)
  582. LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
  583. & + (SHPWRK(2,inode)*xphi1)
  584. LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
  585. & + (SHPWRK(3,inode)*xphi1)
  586. IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  587. & + (SHPWRK(4,inode)*xphi1)
  588. 2522 continue
  589.  
  590. 2521 continue
  591. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  592.  
  593.  
  594. 2520 CONTINUE
  595. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  596.  
  597. c on a construit
  598. C LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  599.  
  600.  
  601. c*********************************************************
  602. c Ajout des fonctions de forme d enrichissement
  603. c + leurs derivees dans repere global
  604. c*********************************************************
  605. CALL SHAPX(XGAU,YGAU,ZGAU,MELE,LV7WRK,NBENRMAX,NBENRJ,
  606. &TLREEL,J,SHPWRK,IRET)
  607.  
  608. c retour a la partie commune aux EF enrichi et non enrichi
  609. 2999 continue
  610.  
  611. C*********************************************************
  612. c write(*,*) 'C CALCUL DE B'
  613. C*********************************************************
  614. call ZERO(BGENE,LHOOK,LRE)
  615. KB=1
  616. C boucle sur tous les Ni
  617. DO 3001 II=1,NBNI
  618.  
  619. BGENE(1,KB) = SHPWRK(2,II)
  620. BGENE(2,KB+1) = SHPWRK(3,II)
  621. BGENE(4,KB) = SHPWRK(3,II)
  622. BGENE(4,KB+1) = SHPWRK(2,II)
  623.  
  624. IF (IDIM.EQ.3) THEN
  625. BGENE(3,KB+2)=SHPWRK(4,II)
  626. BGENE(5,KB) =SHPWRK(4,II)
  627. BGENE(5,KB+2)=SHPWRK(2,II)
  628. BGENE(6,KB+1)=SHPWRK(4,II)
  629. BGENE(6,KB+2)=SHPWRK(3,II)
  630. ENDIF
  631. KB = KB + IDIM
  632. 3001 CONTINUE
  633.  
  634.  
  635. C*********************************************************
  636. c write(*,*) 'C CALCUL DE D (=MATRICE DE HOOKE)'
  637. C*********************************************************
  638. c
  639. c+++++cas DOHOOO
  640. IF (IMAT.EQ.2) THEN
  641. MELVAL=IVAL(1)
  642. IBMN=MIN(J ,IELCHE(/2))
  643. IGMN=MIN(KGAU,IELCHE(/1))
  644. MLREEL=IELCHE(IGMN,IBMN)
  645. SEGACT MLREEL
  646. IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1))
  647. 1 CALL DOHOOO(PROG,LHOOK,DDHOOK)
  648. SEGDES MLREEL
  649.  
  650. c+++++cas DOHMAS
  651. ELSE IF (IMAT.EQ.1) THEN
  652. DO 4001 IM=1,NMATT
  653. IF (IVAL(IM).NE.0) THEN
  654. MELVAL=IVAL(IM)
  655. IBMN=MIN(J ,VELCHE(/2))
  656. IGMN=MIN(KGAU,VELCHE(/1))
  657. if (ibmn.gt.0.and.igmn.gt.0) then
  658. VALMAT(IM)=VELCHE(IGMN,IBMN)
  659. else
  660. VALMAT(IM)=0.D0
  661. endif
  662. ELSE
  663. VALMAT(IM)=0.D0
  664. ENDIF
  665. 4001 CONTINUE
  666. C
  667. CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRET)
  668. IF (IRET.EQ.0) THEN
  669. MOTERR(1:8)=CMATE
  670. MOTERR(9:16)=NOMFR(MFR/2+1)
  671. INTERR(1)=IFOUR
  672. CALL ERREUR(81)
  673. RETURN
  674. ENDIF
  675. ENDIF
  676. C
  677. C*********************************************************
  678. c write(*,*) 'C CALCUL DE Bt.D.B'
  679. C*********************************************************
  680.  
  681. C PRISE EN COMPTE DU POIDS D'INTEGRATION
  682. DJAC = ABS(DJAC) * POIGAU(KGAU)
  683.  
  684. c CAS DES CONTRAINTES PLANES
  685. C recuperation de l'epaisseur: DIM3 donnée facultative du materiau
  686. IF (IFOUR.EQ.-2) THEN
  687. MPTVAL=IVACAR
  688. IF (IVACAR.NE.0) THEN
  689. MELVAL=IVAL(1)
  690. IF (MELVAL.NE.0) THEN
  691. IGMN=MIN(KGAU,VELCHE(/1))
  692. IBMN=MIN(J,VELCHE(/2))
  693. DIM3=VELCHE(IGMN,IBMN)
  694. ELSE
  695. DIM3=1.D0
  696. ENDIF
  697. ENDIF
  698. DJAC=DJAC*DIM3
  699. MPTVAL=IVAMAT
  700. ENDIF
  701.  
  702. CALL ZEROX(REL,NLIGRD,NLIGRP,LRE)
  703.  
  704. CALL BDBST(BGENE,DJAC,DDHOOK,LRE,LHOOK,REL)
  705. * CALL BDBST(BGENE,DJAC,DDHOOK,NLIGRP,LHOOK,REL)
  706. c |-> impossible car BDBST considere qu il sagit de la dimension de bgene
  707. c pour gagner du temps cpu il faudrait qqchose du type:
  708. c CALL BDBSTX(BGENE,DJAC,DDHOOK,NLIGRP,LHOOK,LRE,REL)
  709.  
  710. C*********************************************************
  711. C CUMUL DANS RINT(II,JJ)
  712. C*********************************************************
  713. DO 4201 II=1,NLIGRD
  714. DO 4201 JJ=1,NLIGRP
  715. RINT(II,JJ) = RINT(II,JJ) + REL(II,JJ)
  716. 4201 CONTINUE
  717. c
  718. c
  719. 2500 CONTINUE
  720. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  721. C
  722. C*********************************************************
  723. C write(*,*) 'C STOCKAGE DANS XMATRI de RE(II,JJ)'
  724. c et de XMATRI dans IMATRI
  725. C*********************************************************
  726. C
  727. DO 6001 II=1,NLIGRD
  728.  
  729. DO 6002 JJ=1,NLIGRP
  730.  
  731. c XMATRI.RE(II,JJ,nelrig) = RINT(II,JJ)
  732. c RE(II,JJ,NELRIG) = RINT(II,JJ)
  733. RE(II,JJ,IELRIG) = RINT(II,JJ)
  734.  
  735. c si NON enrichi, on met tout a 0 sauf la diago
  736. c if(mlree1.eq.0) RE(II,JJ) = 0. inutile
  737.  
  738. 6002 CONTINUE
  739.  
  740.  
  741. 6001 CONTINUE
  742. c
  743. c quelques desactivations
  744. do IENR=1,NBENR1
  745. do I=1,NBNN1
  746. mlree1 = TLREEL(IENR,I)
  747. if(mlree1.ne.0) segdes,mlree1
  748. enddo
  749. enddo
  750. c
  751. c
  752. 2001 CONTINUE
  753. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  754.  
  755.  
  756. C
  757. SEGSUP,WRK1,WRK2,MVELCH
  758. SEGSUP,MRACC
  759.  
  760. segdes,MELEME
  761. segdes,MINT1
  762. if(MINT2.ne.0) segdes,MINT2
  763. c
  764. c desactivation des maillages et segment imatri
  765. do ienr=1,(1+NBENR1)
  766. c IPT1 = LOCIRI(1,ienr)
  767. c if(IPT1.ne.0) segdes,IPT1
  768. XMATRI = LOCIRI(4,ienr)
  769. if(XMATRI.ne.0) segdes,XMATRI
  770. enddo
  771.  
  772. c SEGDES dans RIGI1 = IMODEL
  773. C
  774. RETURN
  775. END
  776.  
  777.  
  778.  
  779.  
  780.  

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