Télécharger massx.eso

Retour à la liste

Numérotation des lignes :

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

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