Télécharger sigmax.eso

Retour à la liste

Numérotation des lignes :

  1. C SIGMAX SOURCE BP208322 18/04/12 21:15:13 9802
  2.  
  3. subroutine SIGMAX (MATE,IMAT,NBGMAT,NELMAT,NMATT,CMATE,
  4. & IVAMAT,IMODEL,IREPS2,IVADEP,
  5. & IVASTR,UZDPG,RYDPG,RXDPG,IIPDPG,IRETER)
  6.  
  7. c
  8. C PROCEDURE UTILISEE DANS LE CAS D'ELEMENTS XFEM
  9. c POUR LE CALCUL DE la contrainte (élastique)
  10. C
  11. C
  12. C*********************************************************
  13. C PARTIE DECLARATIVE
  14. C*********************************************************
  15.  
  16. C
  17. IMPLICIT REAL*8(A-H,O-Z)
  18. C
  19. PARAMETER (NDDLMAX=30,NBNIMAX=10)
  20. CTY PARAMETER (NDDLMAX=20,NBNIMAX=10)
  21.  
  22. PARAMETER (NBENRMAX=5)
  23. DIMENSION MLRE(NBENRMAX+1)
  24. C
  25. -INC CCOPTIO
  26. -INC SMCOORD
  27. -INC SMMODEL
  28. -INC SMCHAML
  29. -INC SMINTE
  30. -INC SMELEME
  31. -INC SMLREEL
  32. C
  33. POINTEUR MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE
  34. C
  35. C Segment (type LISTENTI) contenant les informations sur un element
  36. SEGMENT INFO
  37. INTEGER INFELL(JG)
  38. ENDSEGMENT
  39. c
  40. SEGMENT WRK1
  41. REAL*8 XE(3,NBBB)
  42. REAL*8 DDHOOK(LHOOK,LHOOK)
  43. REAL*8 XDDL(LRE),XSTRS(LHOOK)
  44. ENDSEGMENT
  45. c
  46. SEGMENT WRK2
  47. REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE)
  48. c REAL*8 LV7WRK(NBENRMA2,2,6,NBNO)
  49. REAL*8 LV7WRK(NBENRMA2,2,6,NBBB)
  50. ENDSEGMENT
  51. C
  52. SEGMENT,MVELCH
  53. REAL*8 VALMAT(NV1)
  54. ENDSEGMENT
  55. c
  56. SEGMENT MPTVAL
  57. INTEGER IPOS(NS),NSOF(NS)
  58. INTEGER IVAL(NCOSOU)
  59. CHARACTER*16 TYVAL(NCOSOU)
  60. ENDSEGMENT
  61. C
  62. SEGMENT MRACC
  63. INTEGER TLREEL(NBENRMA2,NBI)
  64. ENDSEGMENT
  65. C
  66. SEGMENT MTRACE
  67. REAL*8 TRACE(3,NBPTEL)
  68. ENDSEGMENT
  69. c
  70. CHARACTER*8 CMATE
  71. DIMENSION UDPGE(3)
  72. LOGICAL BDPGE
  73. C
  74. c write (*,*) '#############################'
  75. c write (*,*) '##### DEBUT DE SIGMAX #####'
  76. c write (*,*) '#############################'
  77. C
  78. C*********************************************************
  79. c Introduction du point autour duquel se fait le mouvement
  80. c de la section en defo plane generalisee
  81. C*********************************************************
  82. C En 1D : pas de rotation
  83. BDPGE=.FALSE.
  84. NDPGE=0
  85. XDPGE=0.D0
  86. YDPGE=0.D0
  87. IF (IFOUR.EQ.-3) THEN
  88. BDPGE=.TRUE.
  89. NDPGE=3
  90. UDPGE(1)=UZDPG
  91. UDPGE(2)=RYDPG
  92. UDPGE(3)=RXDPG
  93. SEGACT,MCOORD
  94. IREF=(IIPDPG-1)*(IDIM+1)
  95. XDPGE=XCOOR(IREF+1)
  96. YDPGE=XCOOR(IREF+2)
  97. ELSE IF (IDIM.EQ.1) THEN
  98. IF (IFOUR.EQ. 7 .OR. IFOUR.EQ. 8 .OR. IFOUR.EQ. 9 .OR.
  99. . IFOUR.EQ.10 .OR. IFOUR.EQ.14) THEN
  100. BDPGE=.TRUE.
  101. NDPGE=1
  102. UDPGE(1)=UZDPG
  103. ELSE IF (IFOUR.EQ.11) THEN
  104. BDPGE=.TRUE.
  105. NDPGE=2
  106. UDPGE(1)=UZDPG
  107. UDPGE(2)=RXDPG
  108. ENDIF
  109. ENDIF
  110. C
  111. NV1=NMATT
  112. SEGINI,MVELCH
  113. C
  114. C
  115. C*********************************************************
  116. c
  117. C RECUPERATION + ACTIVATIONS + VALEURS UTILES
  118. c
  119. C*********************************************************
  120. C SEGACT MMODEL,IMODEL deja activé dans epsi1
  121. c
  122. C++++ Activation au cas ou ++++++++++++++++++++++++++++++
  123. SEGACT MCOORD
  124.  
  125. C++++ Recup + Activation de la geometrie ++++++++++++++++
  126. MELEME= IMAMOD
  127. c SEGACT MELEME deja activé dans epsi1
  128.  
  129.  
  130. C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++
  131. c + OBTENUES PAR ELQUOI DANS RIGI1 PENDANT PHASE 1
  132. C segment INFO deja actif dans RIGI1
  133. c + rigi1 n appelle pas elquoi, c'est modeli qui l'a fait
  134. c mais du coup, on na pas de segment minte
  135. c (car depend de si pt de g pour rigi, pour sigma....)
  136. c c'est + simple de rappeler elquoi ici
  137. MELE = NEFMOD
  138. call elquoi(MELE,0,3,IPINF,IMODEL)
  139. INFO = IPINF
  140. c MELE = INFELL(1)
  141. c NBPGA2= INFELL(2)
  142. c NBPGAU= INFELL(3)
  143. c NBPGAU= INFELL(4)
  144. c ICARA = INFELL(5)
  145. NGAU1 = INFELL(6)
  146. c LW = INFELL(7)
  147. LRE = INFELL(9)
  148. LHOOK = INFELL(10)
  149. MINT1 = INFELL(11)
  150. segact,MINT1
  151. MINT2 = INFELL(12)
  152. if(MINT2.ne.0) segact,MINT2
  153. MFR = INFELL(13)
  154. IELE = INFELL(14)
  155. NDDL = INFELL(15)
  156. NSTRS = INFELL(16)
  157. c write(6,*)'-> EPSIX infell',(infell(iou),iou=1,16)
  158.  
  159. c + AUTRES INFOS
  160. C nbre de noeuds par element
  161. NBNN1 = NUM(/1)
  162. C nbre d elements
  163. NBEL1 = NUM(/2)
  164.  
  165. c REM: pour se passer du dimensionnement du nbre d'enrichissement dans
  166. c elquoi et le realiser localement , on pourrait ecrire:
  167. c LRE = NDDLMAX*NBNN1
  168. c NDDL= NDDLMAX
  169.  
  170. C sous decoupage et points de Gauss de l element geometrique de base
  171. CTY if(MELE.eq.263) then
  172. if(MELE.EQ.263.OR.MELE.EQ.264) then
  173. NGAU2 = MINT2.POIGAU(/1)
  174. endif
  175. c write(*,*) 'dim de MINT2=6,',(MINT2.SHPTOT(/2)),(MINT2.SHPTOT(/3))
  176. c write(*,*) 'MINT2',(MINT2.QSIGAU(iou),iou=1,NGAU)
  177. c segdes,MINT2
  178.  
  179.  
  180. c nbre maxi de fonction de forme par noeud (fonction std comprise)
  181. c NBNI = NDDL/IDIM inutile!
  182.  
  183.  
  184. C++++ Recup des infos d enrichissement +++++++++++++++++++
  185. c recup du MCHEX1 d'enrichissement
  186. NOBMO1 = IVAMOD(/1)
  187. if(NOBMO1.ne.0) then
  188. do iobmo1=1,NOBMO1
  189. if((TYMODE(iobmo1)).eq.'MCHAML') then
  190. MCHEX1 = IVAMOD(iobmo1)
  191. segact,MCHEX1
  192. if((MCHEX1.TITCHE).eq.'ENRICHIS') then
  193. MCHAM1 = MCHEX1.ICHAML(1)
  194. segact,MCHAM1
  195. segdes,MCHEX1
  196. goto 1000
  197. endif
  198. segdes,MCHEX1
  199. endif
  200. enddo
  201. write(*,*) 'Le modele est vide (absence d enrichissement)'
  202. * return
  203. else
  204. write(*,*) 'Il n y a pas de MCHEML enrichissement dans le Modele'
  205. * return
  206. endif
  207.  
  208. 1000 continue
  209. c niveau d enrichissement(s) du modele (ddls std u exclus)
  210. c NBENR1= 0 si std, 1 si H seul, 2 si H+F1, 3 si H+F1+F2, etc...
  211. if(NOBMO1.ne.0) then
  212. NBENR1= MCHAM1.IELVAL(/1)
  213. else
  214. NBENR1 = 0
  215. endif
  216. c write(*,*) 'niveau d enrichissement(s) du modele',NBENR1
  217. c
  218. C*********************************************************
  219. C INITIALISATIONS...
  220. C*********************************************************
  221. IRETER = 0
  222. c
  223. c preparation des tables avec:
  224.  
  225. if(NBENR1.ne.0) then
  226. do ienr=1,NBENR1
  227. c -le nombre d'inconnues de chaque sous-zone
  228. c determinee depuis le nombre de fonction de forme
  229. c ienr= 1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10)
  230. nbniJ = 2 + ((ienr-1)*4)
  231. MLRE(1+ienr) = IDIM*NBNN1*nbniJ
  232. enddo
  233. endif
  234. C Tables + longues car 1er indice correspond au fontion de forme std
  235. MLRE(1) = IDIM*NBNN1*1
  236.  
  237. if(NBENR1.lt.(NBENRMAX+1)) then
  238. do ienr=(NBENR1+1),(NBENRMAX)
  239. MLRE(1+ienr) = 0
  240. enddo
  241. endif
  242. c
  243. c ...DU SEGMENT WRK1
  244. NBENRMA2 = NBENRMAX
  245. NBBB = NBNN1
  246. SEGINI,WRK1
  247.  
  248. c ...DU SEGMENT WRK2
  249. c NBNO = NBNI
  250. NBNO = LRE/IDIM
  251. SEGINI,WRK2
  252. C
  253. c ...DU SEGMENT MRACC
  254. NBENRMA2 = NBENRMAX
  255. NBI = NBNN1
  256. segini,MRACC
  257. C
  258. C du nombre d erreur sur le noms de composantes
  259. NBERR1=0
  260.  
  261. C*********************************************************
  262. C
  263. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BOUCLE SUR LES ELEMENTS
  264. C
  265. DO 2000 J=1,NBEL1
  266. c write(*,*) '==================================='
  267. c write(*,*) '========= EF',J,' /NBEL1 ========='
  268. C
  269. C
  270. C*********************************************************
  271. C POUR CHAQUE ELEMENT, ...
  272. C*********************************************************
  273. C
  274. C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT
  275. CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE)
  276. C
  277. c
  278. C++++ NBENRJ = niveau d enrichissement de l element ++++
  279. C =0 si EF std =1 si U+H =2 si U+H+F1 =3 si U+H+F1+F2
  280. NBENRJ=0
  281. if(NBENR1.ne.0) then
  282. do IENR=1,NBENR1
  283. MELVA1 = MCHAM1.IELVAL(IENR)
  284. segact,MELVA1
  285. do I=1,NBNN1
  286. mlree1 = MELVA1.IELCHE(I,J)
  287. C on en profite pour remplir MRACC table de raccourcis pour cet element
  288. TLREEL(IENR,I) = mlree1
  289. if(mlree1.ne.0) then
  290. NBENRJ = max(NBENRJ,IENR)
  291. C et on active la listreel
  292. segact,mlree1
  293. endif
  294. enddo
  295. segdes,MELVA1
  296. enddo
  297. endif
  298. if(NBENRMA2.gt.NBENR1) then
  299. do IENR2=(NBENR1+1),NBENRMA2
  300. do I=1,NBNN1
  301. TLREEL(IENR2,I) = 0
  302. enddo
  303. enddo
  304. endif
  305. C
  306. c
  307. c++++ quelques indices pour dimensionner au plus juste
  308. c nbre total de ddl de l'élément considéré
  309. NLIGRD = MLRE(1+NBENRJ)
  310. NLIGRP = MLRE(1+NBENRJ)
  311. c nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim
  312. NBNI = (MLRE(1+NBENRJ)) / IDIM
  313.  
  314. c write(*,*) 'EF',J,' NBENRJ=',NBENRJ,
  315. c &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme'
  316. c
  317. C++++ ON CALCULE XDDL ++++
  318. MPTVAL = IVADEP
  319. NCOSOU = IVAL(/1)
  320.  
  321. INITOT = 0
  322. C-->> BOUCLE SUR LES niveaux d'ENRICHISSEMENTS (0:U 1:A 2:BCDE1 3:BCDE2)
  323. DO IENR=0,NBENRJ
  324. *nbre de fonction(s) de ce niveau d'enrichissement (=1 si std ou H, =4 pour F1,2,...)
  325. if(IENR .le. 1) then
  326. NBNIENR = 1
  327. else
  328. NBNIENR = 4
  329. endif
  330. C---->> BOUCLE SUR LES fonctions de forme de ce type d'enrichissement
  331. DO INI=1,NBNIENR
  332. INITOT = INITOT + 1
  333. C------>> BOUCLE SUR LA DIMENSION
  334. DO 2220 KDIM=1,IDIM
  335. ICOMP = (INITOT-1)*IDIM + KDIM
  336.  
  337. c --cas ou on n'a pas trouvé assez de composantes--
  338. if(ICOMP.GT.NCOSOU) GOTO 2221
  339.  
  340. MELVAL = IVAL(ICOMP)
  341. c --cas ou on a pas trouvé cette composante dans cette zone du
  342. c chpoint solution devenu mchaml --
  343. if(MELVAL.eq.0)then
  344. c Avait on besoin de cette composante?
  345. c oui, si c'est une composante obligatoire
  346. if(IENR.eq.0) goto 2991
  347. c oui, si l'un des noeuds est enrichi
  348. do I=1,NBNN1
  349. if(TLREEL(IENR,I).ne.0) goto 2991
  350. enddo
  351. c non, si c'est facultatif et qu'on n'est pas enrichi -> on saute
  352. goto 2220
  353. c ->AVERTISSEMENT puis on saute
  354. 2991 NBERR1=NBERR1+1
  355. if(IIMPI.lt.1) goto 2220
  356. c write(IOIMP,*) 'PB OPERATEUR SIGM :'
  357. write(IOIMP,991) ICOMP,IENR,INI,KDIM
  358. 991 format(2X,'ABSENCE DANS LE CHPOINT DEPLACEMENT DE LA ',I3,
  359. $ 'ieme COMPOSANTE (enrichissement',I3,
  360. $ ', fonction',I3,', direction ',I3,
  361. $ ') NECESSAIRE POUR L UN DES NOEUDS SUIVANTS :')
  362. write(IOIMP,*)' noeuds :',(NUM(iou,J),iou=1,NBNN1)
  363. goto 2220
  364. endif
  365.  
  366. C---------->> BOUCLE SUR LES NOEUDS
  367. DO I=1,NBNN1
  368. IQ = ((INITOT-1)*NBNN1*IDIM) + ((I-1)*IDIM) + KDIM
  369. XDDL(IQ) = VELCHE(I,J)
  370. ENDDO
  371.  
  372. 2220 CONTINUE
  373. ENDDO
  374. ENDDO
  375.  
  376. c --cas normal (toutes les composantes souhaitees etaient presentes)--
  377. GOTO 2223
  378.  
  379. c --cas ou on n'a pas trouvé assez de composantes--
  380. 2221 CONTINUE
  381. if (IIMPI.ge.1) then
  382. WRITE(IOIMP,2222) J,NCOSOU,ICOMP
  383. 2222 FORMAT(2X,'ATTENTION : ELEMENT ',I6,
  384. & ' LE CHAMP DE DEPLACEMENT CONTIENT ',I3,' COMPOSANTES',
  385. & ' PAR NOEUD AU LIEU DE ',I3,' ATTENDUES')
  386. endif
  387. NDDL=NCOSOU*NBNN1
  388. NBENRJ=IENR
  389.  
  390. 2223 CONTINUE
  391. c --cas ou on a une ou des erreurs--
  392. IF (NBERR1.gt.0.and.J.eq.NBEL1) THEN
  393. write(IOIMP,*) 'OPERATEUR GRAD : ABSENCE DANS LE CHPOINT ',
  394. & 'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE'
  395. ENDIF
  396.  
  397. c
  398. c rem: il serait probablement interessant au niveau du tems cpu
  399. c d'utiliser moins de pts de Gauss lorsque l element est élastique
  400. c On pourrait par ex. utiliser MINT2 = infell(12) qui contient
  401. c le segment d'integration de l'EF std (QUA4 par ex.)
  402. * if((NBENRJ.eq.0).and.(MINT2.ne.0)) then
  403. * MINTE = MINT2
  404. * NBPGAU= NGAU2
  405. * else
  406. MINTE = MINT1
  407. NBPGAU= NGAU1
  408. * endif
  409. c
  410. c pour les def quadratiques
  411. ISDJC=0
  412. NBPTEL=NBPGAU
  413. IF (IREPS2.EQ.1) SEGINI MTRACE
  414. c
  415. C
  416. C*********************************************************
  417. C
  418. C>>>>>>>>>> BOUCLE SUR LES POINTS DE GAUSS
  419. C
  420. DO 2500 KGAU=1,NBPGAU
  421. C
  422. C*********************************************************
  423. C Initialisation à 0
  424. C*********************************************************
  425.  
  426. c ZERO ne serait-il pas facultatif?
  427. CALL ZERO(SHPWRK,6,NBNO)
  428. C
  429. i6zz = 3
  430. IF (IDIM.EQ.3) i6zz = 4
  431. c do ienr7=1,NBENRMAX
  432. do ienr7=1,NBENRJ
  433. do inod7=1,NBNN1
  434. c do i6=1,6
  435. do i6=1,i6zz
  436. LV7WRK(ienr7,1,i6,inod7) = 0
  437. LV7WRK(ienr7,2,i6,inod7) = 0
  438. enddo
  439. enddo
  440. enddo
  441.  
  442.  
  443. c*********************************************************
  444. c Calcul des fonction de forme std dans repere local
  445. c*********************************************************
  446.  
  447. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  448. c (et donc sur les Ni std)
  449. DO 2510 I=1,NBNN1
  450.  
  451. C++++ Calcul des Ni std
  452. c (rappel: 1:Ni 2:Ni,qsi 3:Ni,eta avec i=1,4)
  453. SHPWRK(1,I) = SHPTOT(1,I,KGAU)
  454. SHPWRK(2,I) = SHPTOT(2,I,KGAU)
  455. SHPWRK(3,I) = SHPTOT(3,I,KGAU)
  456. IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)
  457.  
  458. 2510 CONTINUE
  459. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  460.  
  461.  
  462.  
  463. c*********************************************************
  464. c Passage des fonctions de forme std dans repere global
  465. c*********************************************************
  466.  
  467. C++++ CALCUL DES Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  468. CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  469. c if(J.eq.1.and.KGAU.eq.1)
  470. c &write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1)
  471.  
  472. c*********************************************************
  473. c Si on est pas du tout enrichi on peut sauter une grosse partie
  474. c*********************************************************
  475. if(NBENRJ.eq.0) goto 2999
  476.  
  477. c*********************************************************
  478. c Calcul des level set + leurs derivees dans repere global
  479. c*********************************************************
  480.  
  481. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  482. do 2520 ienr=1,NBENRJ
  483.  
  484. c MELVA1=MCHAM1.IELVAL(IENR)
  485. c segact,MELVA1
  486.  
  487. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  488. DO 2521 I=1,NBNN1
  489.  
  490. C++++ Le I eme noeud est-il ienr-enrichi?
  491. mlree1= TLREEL(IENR,I)
  492. if(mlree1.eq.0) goto 2521
  493.  
  494.  
  495. c Calcul du repere local de fissure(=PSI,PHI)
  496. c (rappel: 1,1:psi 1,2:phi 2,1 psi,x ...etc...)
  497. do 2522 inode=1,NBNN1
  498. c pour le H-enrichissement, on n a pas gardé PSI (inutile)
  499. if(ienr.ne.1) then
  500. c valeur de PSI au inode^ieme noeud
  501. xpsi1 = mlree1.PROG(inode)
  502. c qu on multiplie par la valeur de Ni^std au pt de G considéré
  503. LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
  504. & + (SHPWRK(1,inode)*xpsi1)
  505. LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
  506. & + (SHPWRK(2,inode)*xpsi1)
  507. LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
  508. & + (SHPWRK(3,inode)*xpsi1)
  509. IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  510. & + (SHPWRK(4,inode)*xpsi1)
  511. c valeur de PHI au inode^ieme noeud
  512. xphi1 = mlree1.PROG(NBNN1+inode)
  513. else
  514. xphi1 = mlree1.PROG(inode)
  515. endif
  516. LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
  517. & + (SHPWRK(1,inode)*xphi1)
  518. LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
  519. & + (SHPWRK(2,inode)*xphi1)
  520. LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
  521. & + (SHPWRK(3,inode)*xphi1)
  522. IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  523. & + (SHPWRK(4,inode)*xphi1)
  524. 2522 continue
  525.  
  526. 2521 continue
  527. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  528.  
  529.  
  530. 2520 CONTINUE
  531. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  532.  
  533. c on a construit
  534. C LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  535.  
  536. c*********************************************************
  537. c Ajout des fonctions de forme d enrichissement
  538. c + leurs derivees dans repere global
  539. c*********************************************************
  540. CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET)
  541.  
  542. c retour a la partie commune aux EF enrichi et non enrichi
  543. 2999 continue
  544.  
  545. C*********************************************************
  546. C CALCUL DE B'
  547. C*********************************************************
  548. c ZERO ne serait-il pas facultatif?
  549. c call ZERO(BGENE,LHOOK,NLIGRP)
  550. KB=1
  551. C boucle sur tous les Ni
  552. DO 3001 II=1,NBNI
  553.  
  554. BGENE(1,KB) = SHPWRK(2,II)
  555. BGENE(2,KB+1) = SHPWRK(3,II)
  556. BGENE(4,KB) = SHPWRK(3,II)
  557. BGENE(4,KB+1) = SHPWRK(2,II)
  558.  
  559. IF(IDIM.EQ.3) THEN
  560. BGENE(3,KB+2)=SHPWRK(4,II)
  561. BGENE(5,KB)=SHPWRK(4,II)
  562. BGENE(5,KB+2)=SHPWRK(2,II)
  563. BGENE(6,KB+1)=SHPWRK(4,II)
  564. BGENE(6,KB+2)=SHPWRK(3,II)
  565. ENDIF
  566.  
  567. KB = KB + IDIM
  568.  
  569. 3001 CONTINUE
  570. C
  571. c if(J.eq.5.and.KGAU.eq.1) then
  572. c if(KGAU.eq.1) then
  573. c write(*,*) 'BGENE(1,..)=',(BGENE(1,iou),iou=1,2*NBNI)
  574. c write(*,*) 'BGENE(2,..)=',(BGENE(2,iou),iou=1,2*NBNI)
  575. c write(*,*) 'BGENE(4,..)=',(BGENE(3,iou),iou=1,2*NBNI)
  576. c endif
  577.  
  578. C*********************************************************
  579. C CALCUL DE D
  580. C*********************************************************
  581.  
  582. c on cherche les matrices de Hooke
  583. MPTVAL=IVAMAT
  584. IF (IMAT.EQ.2) THEN
  585. MELVAL=IVAL(1)
  586. IBMN=MIN(J ,IELCHE(/2))
  587. IGMN=MIN(KGAU,IELCHE(/1))
  588. MLREEL=IELCHE(IGMN,IBMN)
  589. SEGACT MLREEL
  590. IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1))
  591. 1 CALL DOHOOO(PROG,LHOOK,DDHOOK)
  592. SEGDES MLREEL
  593.  
  594. ELSE IF (IMAT.EQ.1) THEN
  595. DO 9004 IM=1,NMATT
  596. IF (IVAL(IM).NE.0) THEN
  597. MELVAL=IVAL(IM)
  598. IBMN=MIN(J ,VELCHE(/2))
  599. IGMN=MIN(KGAU,VELCHE(/1))
  600. VALMAT(IM)=VELCHE(IGMN,IBMN)
  601. ELSE
  602. VALMAT(IM)=0.D0
  603. ENDIF
  604. 9004 CONTINUE
  605. c
  606. IF (MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  607. * IF (IGAU.LE.NBGMAT)
  608. * 1 CALL DOHMAO(VALMAT,CMATE,IFOUR,IDIM,TXR,XLOC,XGLOB,D1HO,
  609. * 2 ROTH,DDHOOK,LHOOK,1,IRTD)
  610. write(*,*) 'cas non prevu a ce jour pour les EF xfem'
  611. return
  612. ELSE
  613. IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1))
  614. 1 CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRTD)
  615. ENDIF
  616. ENDIF
  617. c
  618. c
  619. C*********************************************************
  620. C CALCUL DE sigma = D*B*u
  621. C*********************************************************
  622. c
  623. CALL DBST(BGENE,DDHOOK,XDDL,(NBNI*IDIM),LHOOK,XSTRS)
  624. c
  625. c calcul des eps 2 (termes quadratiques de la deformation)
  626. c
  627. IF (IREPS2.EQ.1)
  628. 1 CALL DBST2(SHPWRK,DDHOOK,XDDL,XE,NBNI,IFOUR,LHOOK,XSTRS,
  629. 2 TRACE,KGAU,XDPGE,YDPGE,UDPGE,NIFOUR)
  630. c
  631. c remplissage du segment contenant les contraintes
  632. c
  633. MPTVAL=IVASTR
  634. DO 7004 ICOMP=1,LHOOK
  635. MELVAL=IVAL(ICOMP)
  636. IBMN=MIN(J,VELCHE(/2))
  637. VELCHE(KGAU,IBMN)=XSTRS(ICOMP)
  638. 7004 CONTINUE
  639. c
  640. * if(J.eq.5 .and. KGAU.eq.1) then
  641. c if(KGAU.eq.1) then
  642. c write(*,*) J,KGAU,'sigma(..)=',(XSTRS(iou),iou=1,LHOOK)
  643. c endif
  644. c
  645. 2500 CONTINUE
  646. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  647. C
  648. c quelques suppressions
  649. c (ici element non-incompressible=> pas besoin de MTRACE (cf epsi2)
  650. IF (IREPS2.EQ.1) THEN
  651. SEGSUP MTRACE
  652. ENDIF
  653. c
  654. c quelques desactivations
  655. do IENR=1,NBENR1
  656. do I=1,NBNN1
  657. mlree1 = TLREEL(IENR,I)
  658. if(mlree1.ne.0) segdes,mlree1
  659. enddo
  660. enddo
  661. c
  662. 2000 CONTINUE
  663. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  664. c
  665. C*********************************************************
  666. C SUPPRESSION ET DESACTIVATION DE SEGMENTS
  667. C*********************************************************
  668. C
  669. SEGSUP,WRK1,WRK2,MVELCH
  670.  
  671. SEGSUP,MRACC
  672.  
  673. segdes,MINT1
  674. if(MINT2.ne.0) segdes,MINT2
  675. C
  676. RETURN
  677. END
  678.  
  679.  
  680.  
  681.  
  682.  

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