Télécharger rigix.eso

Retour à la liste

Numérotation des lignes :

  1. C RIGIX SOURCE GG250959 17/09/21 21:15:19 9560
  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 Segment contenant l'info DDL a mettre à 0
  95. SEGMENT TBLOQ
  96. INTEGER MBLOQ(3,NBPTB)
  97. INTEGER NBLOQ(3)
  98. ENDSEGMENT
  99. C
  100. C
  101. C*********************************************************
  102. c
  103. C RECUPERATION + ACTIVATIONS + VALEURS UTILES
  104. c
  105. C*********************************************************
  106. C sont deja actifs dans rigi1.eso :
  107. C SEGACT MMODEL,IMODEL,MELVAL
  108. c IVAMAT
  109. MPTVAL=IVAMAT
  110.  
  111. c
  112. C++++ Activation au cas ou ++++++++++++++++++++++++++++++
  113. SEGACT MCOORD
  114.  
  115. C++++ Recup + Activation de la geometrie ++++++++++++++++
  116. MELEME= IMAMOD
  117. SEGACT MELEME
  118.  
  119. C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++
  120. c + OBTENUES PAR ELQUOI DANS RIGI1 PENDANT PHASE 1
  121. C segment INFO deja actif dans RIGI1
  122. c + rigi1 n appelle pas elquoi, c'est modeli qui l'a fait
  123. c mais du coup, on na pas de segment minte
  124. c (car depend de si pt de g pour rigi, pour sigma....)
  125. c c'est + simple de rappeler elquoi ici
  126. MELE = NEFMOD
  127. call elquoi(MELE,0,3,IPINF,IMODEL)
  128. INFO = IPINF
  129. c MELE = INFELL(1)
  130. c NBPGA2= INFELL(2)
  131. c NBPGAU= INFELL(3)
  132. c NBPGAU= INFELL(4)
  133. c ICARA = INFELL(5)
  134. NGAU1 = INFELL(6)
  135. c LW = INFELL(7)
  136. LRE = INFELL(9)
  137. LHOOK = INFELL(10)
  138. MINT1 = INFELL(11)
  139. segact,MINT1
  140. MINT2 = INFELL(12)
  141. c if(MINT2.ne.0) segact,MINT2
  142. MFR = INFELL(13)
  143. IELE = INFELL(14)
  144. NDDL = INFELL(15)
  145. NSTRS = INFELL(16)
  146. c write(ioimp,*)'-> RIGIX infell',(infell(iou),iou=1,16)
  147.  
  148. c + AUTRES INFOS
  149. C nbre de noeuds par element
  150. NBNN1 = NUM(/1)
  151. C nbre d elements
  152. NBEL1 = NUM(/2)
  153.  
  154. c REM: pour se passer du dimensionnement du nbre d'enrichissement dans
  155. c elquoi et le realiser localement , on pourrait ecrire:
  156. c LRE = NDDLMAX*NBNN1
  157. c NDDL= NDDLMAX
  158.  
  159. C sous decoupage et points de Gauss de l element geometrique de base
  160. CYT if(MELE.eq.263) then
  161. cbp if(MELE.eq.263.or.MELE.eq.264) then
  162. if (MINT2.ne.0) then
  163. segact,MINT2
  164. NGAU2 = MINT2.POIGAU(/1)
  165. else
  166. NGAU2=0
  167. endif
  168. c write(*,*) 'dim de MINT2=6,',(MINT2.SHPTOT(/2)),(MINT2.SHPTOT(/3))
  169. c write(*,*) 'MINT2',(MINT2.QSIGAU(iou),iou=1,NGAU2)
  170. c segdes,MINT2
  171.  
  172.  
  173. c nbre maxi de fonction de forme par noeud (fonction std comprise)
  174. c NBNI = NDDL/IDIM inutile!
  175.  
  176.  
  177. C++++ Recup des infos d enrichissement +++++++++++++++++++
  178. c recup du MCHEX1 d'enrichissement
  179. NBENR1=0
  180. MCHAM1=0
  181. NOBMO1=IVAMOD(/1)
  182. if (NOBMO1.ne.0) then
  183. do iobmo1=1,NOBMO1
  184. if ((TYMODE(iobmo1)).eq.'MCHAML') then
  185. MCHEX1 = IVAMOD(iobmo1)
  186. segact,MCHEX1
  187. if ((MCHEX1.TITCHE).eq.'ENRICHIS') then
  188. MCHAM1 = MCHEX1.ICHAML(1)
  189. segact,MCHAM1
  190. segdes,MCHEX1
  191. goto 1000
  192. endif
  193. segdes,MCHEX1
  194. endif
  195. enddo
  196. write(ioimp,*) 'Le modele est vide (absence d enrichissement)'
  197. * return
  198. else
  199. write(ioimp,*) 'Aucun MCHAML enrichissement dans le Modele'
  200. * return
  201. endif
  202.  
  203. 1000 continue
  204. c niveau d enrichissement(s) du modele (ddls std u exclus)
  205. c NBENR1= 0 si std, 1 si H seul, 2 si H+F1, 3 si H+F1+F2, etc...
  206. if (MCHAM1.ne.0) NBENR1=MCHAM1.IELVAL(/1)
  207. c* write(ioimp,*) 'niveau d enrichissement(s) du modele',NBENR1
  208.  
  209. C*********************************************************
  210. c
  211. c PREPARATION DES OBJETS RESULTATS
  212. c
  213. C*********************************************************
  214. C
  215. C++++ RECHERCHE DES NOMS D'INCONNUES ET DES DUAUX
  216. c MOTINC et MODUA en dur pour l instant
  217.  
  218. C++++ REMPLISSAGE DE DESCR de taille maxi
  219. c (on ajustera en enlevant si besoin)
  220. NLIGRP = LRE
  221. NLIGRD = LRE
  222. SEGINI DESCR
  223. IPDSCR = DESCR
  224. C
  225. KINC = 0
  226. C----->> BOUCLE SUR LES fonctions de Forme
  227. DO IENR=1,NBNIMAX
  228. C------->> BOUCLE SUR LES NOEUDS
  229. DO I=1,NBNN1
  230. C--------->> BOUCLE SUR LA DIMENSION
  231. DO KDIM=1,IDIM
  232. KINC = KINC + 1
  233. CYT LISINC(KINC) = MOTINC((IDIM*(IENR-1))+KDIM)
  234. CYT LISDUA(KINC) = MOTDUA((IDIM*(IENR-1))+KDIM)
  235. LISINC(KINC) = MOTINC((3*(IENR-1))+KDIM)
  236. LISDUA(KINC) = MOTDUA((3*(IENR-1))+KDIM)
  237. NOELEP(KINC) = I
  238. NOELED(KINC) = I
  239. ENDDO
  240. ENDDO
  241. ENDDO
  242.  
  243. IF (KINC.NE.LRE) THEN
  244. WRITE(ioimp,*) 'IL Y A UNE ERREUR DANS DESCR'
  245. WRITE(ioimp,*) 'KINC , LRE = ',KINC,LRE
  246. RETURN
  247. ENDIF
  248. C
  249. SEGDES DESCR
  250. C
  251. C
  252. C++++ INITIALISATIONS...
  253. C
  254. c ... des tables LOCIRI et MLRE
  255. c MLRE contient le nombre d'inconnues de chaque sous-zone
  256. c determinee depuis le nombre de fonctions de forme
  257. c ienr= 1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10)
  258. if (NBENR1.ne.0) then
  259. do ienr=1,NBENR1
  260. nbniJ = 2 + ((ienr-1)*4)
  261. MLRE(1+ienr) = IDIM*NBNN1*nbniJ
  262. c -LOCIRI
  263. LOCIRI(5,1+ienr)= NIFOUR
  264. enddo
  265. endif
  266. C Tables + longues car 1er indice correspond au fontion de forme std
  267. MLRE(1) = IDIM*NBNN1*1
  268. LOCIRI(5,1)= NIFOUR
  269. c on complete avec des 0
  270. if (NBENR1.lt.(NBENRMAX+1)) then
  271. do ienr=(NBENR1+1),(NBENRMAX)
  272. MLRE(1+ienr) = 0
  273. enddo
  274. endif
  275. c
  276. c ...DU SEGMENT WRK1
  277. NBENRMA2 = NBENRMAX
  278. NBBB = NBNN1
  279. SEGINI,WRK1
  280.  
  281. c ...DU SEGMENT WRK2
  282. c NBNO = NBNI
  283. NBNO = LRE/IDIM
  284. SEGINI,WRK2
  285.  
  286. c ...DU SEGMENT MVELCH
  287. NV1 = NMATT
  288. SEGINI,MVELCH
  289. C
  290. c ...DU SEGMENT MRACC
  291. NBENRMA2 = NBENRMAX
  292. NBI = NBNN1
  293. NBELEM = NUM(/2)
  294. segini,MRACC
  295. C
  296. C
  297.  
  298. C*********************************************************
  299. C
  300. C>>>>>>>>>>>>>>>>>>>>>>>>>>> 1ere BOUCLE SUR LES ELEMENTS
  301. C
  302. NBENR = 0
  303. DO 2000 J=1,NBEL1
  304. c write(ioimp,*) '===boucle1=== EF',J,' /NBEL1 ================'
  305. C
  306. C
  307. C*********************************************************
  308. C POUR CHAQUE ELEMENT, ...
  309. C*********************************************************
  310. C
  311. C++++ NBENRJ = niveau d enrichissement de l element ++++
  312. C =0 si EF std =1 si U+H =2 si U+H+F1 =3 si U+H+F1+F2
  313. NBENRJ=0
  314. if (NBENR1.ne.0) then
  315. do IENR=1,NBENR1
  316. MELVA1 = MCHAM1.IELVAL(IENR)
  317. segact,MELVA1
  318. do I=1,NBNN1
  319. mlree1 = MELVA1.IELCHE(I,J)
  320. if(mlree1.ne.0) NBENRJ=max(NBENRJ,IENR)
  321. enddo
  322. segdes,MELVA1
  323. enddo
  324. endif
  325. NBENR=max(NBENRJ,NBENR)
  326. c
  327. c++++ Selon le niveau d enrichissement NBENRJ...
  328. C
  329. c++++ ...on ajuste la 3eme dimension de XMATRI.RE
  330. c (=nbre d element de ce niveau d enrichissement)
  331. C on la stocke temporairement dans LOCIRI(4, ) meme si c est pas propre
  332. LOCIRI(4,1+NBENRJ) = LOCIRI(4,1+NBENRJ) + 1
  333. c
  334. C++++ ...on copie cet element a la geo correspondante
  335. IPT1 = LOCIRI(1,1+NBENRJ)
  336. c si elle n existe pas, il faut la creer
  337. if (IPT1.eq.0) then
  338. NBNN = NBNN1
  339. NBELEM = 1
  340. NBSOUS = 0
  341. NBREF = 0
  342. segini,IPT1
  343. IPT1.ITYPEL = ITYPEL
  344. LOCIRI(1,1+NBENRJ)=IPT1
  345. c si elle existe, il faut l agrandir
  346. else
  347. NBNN = NBNN1
  348. NBELEM = (IPT1.NUM(/2)) + 1
  349. NBSOUS = 0
  350. NBREF = 0
  351. segadj,IPT1
  352. endif
  353. c copie en cours
  354. do I=1,NBNN1
  355. IPT1.NUM(I,NBELEM) = NUM(I,J)
  356. enddo
  357.  
  358. c on retient la position NBELEM de l element J
  359. MELRIG(J)= NBELEM
  360.  
  361. C++++ ...on crée le descr qui va bien si il n existe pas
  362. DES1 = LOCIRI(3,1+NBENRJ)
  363. if (DES1.eq.0) then
  364. segini,DES1=DESCR
  365. NLIGRP = MLRE(1+NBENRJ)
  366. NLIGRD = MLRE(1+NBENRJ)
  367. segadj,DES1
  368. LOCIRI(3,1+NBENRJ) = DES1
  369. segdes,DES1
  370. endif
  371.  
  372. 2000 CONTINUE
  373. C FIN DE 1ere BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  374. C
  375. C*********************************************************
  376.  
  377.  
  378. C*********************************************************
  379. C POUR CHAQUE NIVEAU D ENRICHISSEMENT, ...
  380. C*********************************************************
  381. do IENR=1,(NBENR1+1)
  382. c attention : IENR = IENRvrai + 1
  383. NELRIG = LOCIRI(4,IENR)
  384. C++++ ...ON CRÉE XMATRI DE LA BONNE TAILLE TOUT DE SUITE
  385. if (NELRIG.ne.0) then
  386. NLIGRP = MLRE(IENR)
  387. NLIGRD = MLRE(IENR)
  388. segini,XMATRI
  389. LOCIRI(4,IENR)=XMATRI
  390. C++++ ...ON DESACTIVE IPT1
  391. IPT1 = LOCIRI(1,IENR)
  392. segdes,IPT1
  393. endif
  394. enddo
  395.  
  396.  
  397. C initialisation du tableau des ddl a mettre à zéro
  398. NBPTB= (XCOOR(/1))/(IDIM+1)
  399. SEGINI TBLOQ
  400.  
  401. C++++ TBLOQ.MBLOQ(NBENRJ, INUM) = faut il mettre a 0 les ddl ++++
  402. C =0 si pas encore passé dans les ddl
  403. C =1 si déja passé dans les ddl pas de mise à 0
  404. c =2 si déja passé dans les ddl mise à 0 nécéssaire
  405. C ++++ TBLOQ.NBLOQ(NBENRJ) Compteur de ddl de type NBENRJ à mettre à 0
  406.  
  407.  
  408. C*********************************************************
  409. C
  410. C>>>>>>>>>>>>>>>>>>>>>>>>>>> 2eme BOUCLE SUR LES ELEMENTS
  411. C
  412. DO 2001 J=1,NBEL1
  413. c write(*,*) '===boucle 2=== EF',J,' /NBEL1 ================'
  414. C
  415. C*********************************************************
  416. C POUR CHAQUE ELEMENT, ...
  417. C*********************************************************
  418. C
  419. C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT
  420. CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE)
  421.  
  422. C++++ NBENRJ = niveau d enrichissement de l element ++++
  423. C =0 si EF std =1 si U+H =2 si U+H+F1 =3 si U+H+F1+F2
  424. NBENRJ=0
  425. if (NBENR1.ne.0) then
  426. do IENR=1,(NBENR1)
  427. MELVA1 = MCHAM1.IELVAL(IENR)
  428. segact,MELVA1
  429. do I=1,NBNN1
  430. mlree1 = MELVA1.IELCHE(I,J)
  431. C on en profite pour remplir MRACC table de raccourcis pour cet element
  432. TLREEL(IENR,I) = mlree1
  433. if (mlree1.ne.0) then
  434. NBENRJ = max(NBENRJ,IENR)
  435. C et on active la listreel
  436. segact,mlree1
  437. endif
  438. enddo
  439. segdes,MELVA1
  440. enddo
  441. endif
  442. if (NBENRMA2.gt.NBENR1) then
  443. do IENR2=(NBENR1+1),NBENRMA2
  444. do I=1,NBNN1
  445. TLREEL(IENR2,I) = 0
  446. enddo
  447. enddo
  448. endif
  449. C
  450. c++++ on recupere XMATRI
  451. XMATRI = LOCIRI(4,1+NBENRJ)
  452.  
  453. c++++ quelques indices pour dimensionner au plus juste
  454. c nbre total de ddl de l'élément considéré
  455. NLIGRD = MLRE(1+NBENRJ)
  456. NLIGRP = MLRE(1+NBENRJ)
  457. c nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim
  458. NBNI = (MLRE(1+NBENRJ)) / IDIM
  459. c position de cet element pour cet enrichissement
  460. IELRIG = MELRIG(J)
  461. c write(*,*) 'EF',J,' NBENRJ=',NBENRJ,
  462. c &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme'
  463. c
  464.  
  465. C++++ ON MET A ZERO LA SOMME QUE L ON VA CALCULER
  466. c CALL ZERO(RINT,LRE,NLIGRP)
  467. * CALL ZERO(RINT,NLIGRD,NLIGRP)
  468. c |-> impossible car ZERO considere qu il sagit de la dimension de SHPWRK
  469. CALL ZEROX(RINT,NLIGRD,NLIGRP,LRE)
  470. c
  471. c REM : il est interessant au niveau du tems cpu d'utiliser moins de pts
  472. c de Gauss lorsque l element n est pas enrichi car cela n'apporte
  473. c rien de plus.
  474. c On utilise MINT2 = infell(12) qui contient le segment d'integration
  475. C de l'EF std (QUA4 par ex. pour XQ4R)
  476. if ((NBENRJ.eq.0).and.(MINT2.ne.0)) then
  477. MINTE = MINT2
  478. NBPGAU= NGAU2
  479. else
  480. MINTE = MINT1
  481. NBPGAU= NGAU1
  482. endif
  483.  
  484. C*********************************************************
  485. C
  486. C>>>>>>>>>> BOUCLE SUR LES POINTS DE GAUSS
  487. C
  488. DO 2500 KGAU=1,NBPGAU
  489. c
  490. c write(ioimp,*) '-pt de G ',KGAU,' / ',NBPGAU,'----'
  491. C
  492. C*********************************************************
  493. C Initialisation à 0
  494. C*********************************************************
  495.  
  496. c CALL ZERO(SHPWRK,6,NBNO)
  497. CALL ZERO(SHPWRK,6,NBNI)
  498.  
  499. XGAU = 0.D0
  500. YGAU = 0.D0
  501. ZGAU = 0.D0
  502. C
  503. i6zz = 3
  504. IF (IDIM.EQ.3) i6zz = 4
  505. C
  506. c do ienr7=1,NBENRMAX
  507. do ienr7=1,NBENRJ
  508. do inod7=1,NBNN1
  509. c do i6=1,6
  510. CYT do i6=1,3
  511. do i6=1,i6zz
  512. LV7WRK(ienr7,1,i6,inod7) = 0.D0
  513. LV7WRK(ienr7,2,i6,inod7) = 0.D0
  514. enddo
  515. enddo
  516. enddo
  517.  
  518.  
  519. c*********************************************************
  520. c Calcul des fonction de forme std dans repere local
  521. c*********************************************************
  522.  
  523. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  524. c (et donc sur les Ni std)
  525. DO 2510 I=1,NBNN1
  526.  
  527. C++++ Calcul des Ni std
  528. c (rappel: 1:Ni 2:Ni,qsi 3:Ni,eta avec i=1,4)
  529. SHPWRK(1,I) = SHPTOT(1,I,KGAU)
  530. SHPWRK(2,I) = SHPTOT(2,I,KGAU)
  531. SHPWRK(3,I) = SHPTOT(3,I,KGAU)
  532. IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)
  533.  
  534. C++++ Calcul des coordonnees dans le repere global
  535. XGAU = XGAU + (SHPWRK(1,I)*XE(1,I))
  536. YGAU = YGAU + (SHPWRK(1,I)*XE(2,I))
  537. IF (IDIM.EQ.3) ZGAU = ZGAU + (SHPWRK(1,I)*XE(3,I))
  538.  
  539. 2510 CONTINUE
  540. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  541.  
  542.  
  543.  
  544. c*********************************************************
  545. c Passage des fonctions de forme std dans repere global
  546. c*********************************************************
  547.  
  548. C++++ CALCUL DES Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  549. CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  550. c if(J.eq.1.and.KGAU.eq.1)
  551. c &write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1)
  552.  
  553. c*********************************************************
  554. c Si on est pas du tout enrichi on peut sauter une grosse partie
  555. c*********************************************************
  556. if(NBENRJ.eq.0) goto 2999
  557.  
  558. c*********************************************************
  559. c Calcul des level set + leurs derivees dans repere global
  560. c*********************************************************
  561.  
  562. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  563. do 2520 ienr=1,NBENRJ
  564.  
  565. c MELVA1=MCHAM1.IELVAL(IENR)
  566. c segact,MELVA1
  567.  
  568. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  569. DO 2521 I=1,NBNN1
  570.  
  571. C++++ Le I eme noeud est-il ienr-enrichi?
  572. mlree1= TLREEL(IENR,I)
  573. if(mlree1.eq.0) goto 2521
  574.  
  575.  
  576. c Calcul du repere local de fissure(=PSI,PHI)
  577. c (rappel: 1,1:psi 1,2:phi 2,1 psi,x ...etc...)
  578. do 2522 inode=1,NBNN1
  579. c pour le H-enrichissement, on n a pas gardé PSI (inutile)
  580. if (ienr.ne.1) then
  581. c valeur de PSI au inode^ieme noeud
  582. xpsi1 = mlree1.PROG(inode)
  583. c qu on multiplie par la valeur de Ni^std au pt de G considéré
  584. LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
  585. & + (SHPWRK(1,inode)*xpsi1)
  586. LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
  587. & + (SHPWRK(2,inode)*xpsi1)
  588. LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
  589. & + (SHPWRK(3,inode)*xpsi1)
  590. IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  591. & + (SHPWRK(4,inode)*xpsi1)
  592. c valeur de PHI au inode^ieme noeud
  593. xphi1 = mlree1.PROG(NBNN1+inode)
  594. else
  595. xphi1 = mlree1.PROG(inode)
  596. endif
  597. LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
  598. & + (SHPWRK(1,inode)*xphi1)
  599. LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
  600. & + (SHPWRK(2,inode)*xphi1)
  601. LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
  602. & + (SHPWRK(3,inode)*xphi1)
  603. IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  604. & + (SHPWRK(4,inode)*xphi1)
  605. 2522 continue
  606.  
  607. 2521 continue
  608. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  609.  
  610.  
  611. 2520 CONTINUE
  612. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  613.  
  614. c on a construit
  615. C LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  616.  
  617.  
  618. c*********************************************************
  619. c Ajout des fonctions de forme d enrichissement
  620. c + leurs derivees dans repere global
  621. c*********************************************************
  622.  
  623. CALL SHAPX(XGAU,YGAU,ZGAU,MELE,LV7WRK,NBENRMAX,NBENRJ,
  624. &TLREEL,J,SHPWRK,IRET)
  625.  
  626. c retour a la partie commune aux EF enrichi et non enrichi
  627. 2999 continue
  628.  
  629. C*********************************************************
  630. c write(*,*) 'C CALCUL DE B'
  631. C*********************************************************
  632. call ZERO(BGENE,LHOOK,LRE)
  633. KB=1
  634. C boucle sur tous les Ni
  635. DO 3001 II=1,NBNI
  636.  
  637. BGENE(1,KB) = SHPWRK(2,II)
  638. BGENE(2,KB+1) = SHPWRK(3,II)
  639. BGENE(4,KB) = SHPWRK(3,II)
  640. BGENE(4,KB+1) = SHPWRK(2,II)
  641.  
  642. IF (IDIM.EQ.3) THEN
  643. BGENE(3,KB+2)=SHPWRK(4,II)
  644. BGENE(5,KB) =SHPWRK(4,II)
  645. BGENE(5,KB+2)=SHPWRK(2,II)
  646. BGENE(6,KB+1)=SHPWRK(4,II)
  647. BGENE(6,KB+2)=SHPWRK(3,II)
  648. ENDIF
  649. KB = KB + IDIM
  650. 3001 CONTINUE
  651.  
  652.  
  653. C*********************************************************
  654. c write(*,*) 'C CALCUL DE D (=MATRICE DE HOOKE)'
  655. C*********************************************************
  656. c
  657. c+++++cas DOHOOO
  658. IF (IMAT.EQ.2) THEN
  659. MELVAL=IVAL(1)
  660. IBMN=MIN(J ,IELCHE(/2))
  661. IGMN=MIN(KGAU,IELCHE(/1))
  662. MLREEL=IELCHE(IGMN,IBMN)
  663. SEGACT MLREEL
  664. IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1))
  665. 1 CALL DOHOOO(PROG,LHOOK,DDHOOK)
  666. SEGDES MLREEL
  667.  
  668. c+++++cas DOHMAS
  669. ELSE IF (IMAT.EQ.1) THEN
  670. DO 4001 IM=1,NMATT
  671. IF (IVAL(IM).NE.0) THEN
  672. MELVAL=IVAL(IM)
  673. IBMN=MIN(J ,VELCHE(/2))
  674. IGMN=MIN(KGAU,VELCHE(/1))
  675. if (ibmn.gt.0.and.igmn.gt.0) then
  676. VALMAT(IM)=VELCHE(IGMN,IBMN)
  677. else
  678. VALMAT(IM)=0.D0
  679. endif
  680. ELSE
  681. VALMAT(IM)=0.D0
  682. ENDIF
  683. 4001 CONTINUE
  684. C
  685. CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRET)
  686. IF (IRET.EQ.0) THEN
  687. MOTERR(1:8)=CMATE
  688. MOTERR(9:16)=NOMFR(MFR/2+1)
  689. INTERR(1)=IFOUR
  690. CALL ERREUR(81)
  691. RETURN
  692. ENDIF
  693. ENDIF
  694. C
  695. C*********************************************************
  696. c write(*,*) 'C CALCUL DE Bt.D.B'
  697. C*********************************************************
  698.  
  699. C PRISE EN COMPTE DU POIDS D'INTEGRATION
  700. DJAC = ABS(DJAC) * POIGAU(KGAU)
  701.  
  702. c CAS DES CONTRAINTES PLANES
  703. C recuperation de l'epaisseur: DIM3 donnée facultative du materiau
  704. IF (IFOUR.EQ.-2) THEN
  705. MPTVAL=IVACAR
  706. IF (IVACAR.NE.0) THEN
  707. MELVAL=IVAL(1)
  708. IF (MELVAL.NE.0) THEN
  709. IGMN=MIN(KGAU,VELCHE(/1))
  710. IBMN=MIN(J,VELCHE(/2))
  711. DIM3=VELCHE(IGMN,IBMN)
  712. ELSE
  713. DIM3=1.D0
  714. ENDIF
  715. ENDIF
  716. DJAC=DJAC*DIM3
  717. MPTVAL=IVAMAT
  718. ENDIF
  719.  
  720. CALL ZEROX(REL,NLIGRD,NLIGRP,LRE)
  721.  
  722. CALL BDBST(BGENE,DJAC,DDHOOK,LRE,LHOOK,REL)
  723. * CALL BDBST(BGENE,DJAC,DDHOOK,NLIGRP,LHOOK,REL)
  724. c |-> impossible car BDBST considere qu il sagit de la dimension de bgene
  725. c pour gagner du temps cpu il faudrait qqchose du type:
  726. c CALL BDBSTX(BGENE,DJAC,DDHOOK,NLIGRP,LHOOK,LRE,REL)
  727.  
  728. C*********************************************************
  729. C CUMUL DANS RINT(II,JJ)
  730. C*********************************************************
  731. DO 4201 II=1,NLIGRD
  732. DO 4201 JJ=1,NLIGRP
  733. RINT(II,JJ) = RINT(II,JJ) + REL(II,JJ)
  734. 4201 CONTINUE
  735. c
  736. c
  737. 2500 CONTINUE
  738. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  739. C
  740. C*********************************************************
  741. C write(*,*) 'C STOCKAGE DANS XMATRI de RE(II,JJ)'
  742. c et de XMATRI dans IMATRI
  743. C*********************************************************
  744. C
  745. DO 6001 II=1,NLIGRD
  746.  
  747. DO 6002 JJ=1,NLIGRP
  748.  
  749. c XMATRI.RE(II,JJ,nelrig) = RINT(II,JJ)
  750. c RE(II,JJ,NELRIG) = RINT(II,JJ)
  751. RE(II,JJ,IELRIG) = RINT(II,JJ)
  752.  
  753. c si NON enrichi, on met tout a 0 sauf la diago
  754. c if(mlree1.eq.0) RE(II,JJ) = 0. inutile
  755.  
  756. 6002 CONTINUE
  757.  
  758.  
  759.  
  760. C**********************************************************************
  761. C On cherche les noeuds qui ne sont pas effectivement enrichis
  762. Cpour forcer à 0 les DDL correspondants
  763. C**********************************************************************
  764.  
  765. c on trouve le type de fonction de ce ddl: IENR=0 si U, =1 si A,
  766. c =2 si B1, =3 si C1, = 4 si D1, =5 si E1, =6 si B2, ... =9 si E2
  767. IENR = ((II-1)/IDIM) / NBNN1
  768. c on trouve le niveau d enrichissement KENR de ce ddl si NONstd
  769. if (IENR.eq.0) then
  770. go to 6011
  771. elseif(IENR.ge.2) then
  772. KENR = ((IENR - 2) / 4) + 2
  773. c ci dessus: 4 represente le nombre de fonction de la base d'enrichissement
  774. c et 2 est le decalage du a U et H
  775. else
  776. KENR = IENR
  777. endif
  778. c on trouve le noeud correspondant au ddl
  779. CYT INODE = ((II+1)/IDIM) - ((IENR)*NBNN1)
  780. INODE = 1 + ((II-1)/IDIM) - ((IENR)*NBNN1)
  781. c numero global du noeud
  782. INUM = NUM(INODE , J)
  783. c est ont déja passé dans ce ddl ?
  784. c write(*,*) 'INUM', INUM, 'Kenr', KENR
  785. c write(*,*) 'mlree1', mlree1,'Mbloq',Tbloq.Mbloq(KENR, INUM)
  786. if (Tbloq.Mbloq(KENR, INUM).gt.0) GOTO 6011
  787. c est-ce un noeud vraiment enrichi?
  788. mlree1 = TLREEL(KENR,INODE)
  789. Tbloq.Mbloq(KENR, INUM)=Tbloq.Mbloq(KENR, INUM)+1
  790. if (mlree1.eq.0) then
  791. Tbloq.Mbloq(KENR, INUM)=Tbloq.Mbloq(KENR, INUM)+1
  792. Tbloq.Nbloq(KENR) = Tbloq.Nbloq(KENR)+1
  793. c XMATRI.RE(II,II,nelrig) = 1.D0
  794. c RE(II,II,NELRIG) = 1.D0
  795. c RE(II,II,IELRIG) = 1.D0
  796. endif
  797. c write(*,*) 'Nouveau Mbloq', Tbloq.Mbloq(KENR, INUM)
  798.  
  799.  
  800. 6011 CONTINUE
  801. 6001 CONTINUE
  802. c
  803. c quelques desactivations
  804. do IENR=1,NBENR1
  805. do I=1,NBNN1
  806. mlree1 = TLREEL(IENR,I)
  807. if(mlree1.ne.0) segdes,mlree1
  808. enddo
  809. enddo
  810. c
  811. c
  812. 2001 CONTINUE
  813. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  814.  
  815. C*********************************************************
  816. C Creation des matricres de bloquage
  817. C pour les DDL non effectivement enrichis
  818. C*********************************************************
  819.  
  820. NLIGRD = 2
  821. NLIGRP = 2
  822.  
  823. NBNN = 2
  824. NBPTS1=NBPTB
  825.  
  826. C Boucle sur les types d'enrichissement <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  827. NBRIG = 1 + NBENR
  828. IDDL = 3
  829. DO 2002 IENR =1, NBENR1
  830. C Maillage des noeuds à bloquer
  831. C nombre de noeuds a bloquer
  832. if (IENR.eq.1) then
  833. NBELEM = TBLOQ.NBLOQ(1)
  834. elseif (IENR.EQ.2) then
  835. NBELEM = TBLOQ.NBLOQ(2)
  836. else
  837. NBELEM = TBLOQ.NBLOQ(3)
  838. endif
  839. NFON=1
  840. IF (IENR.gt.1) NFON=4
  841. C Boucle sur les fonctions de formes d'enrichissement
  842. DO 2032 IFON = 1 , NFON
  843. C Boucle sur les composantes <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  844. DO 2022 ICOMP = 1, IDIM
  845.  
  846. NELRIG = NBELEM
  847. IDDL=IDDL+1
  848. C si aucun noeud a bloquer on saute cette composante
  849. if (NELRIG.EQ.0) goto 2002
  850. NBRIG = NBRIG + 1
  851.  
  852. SEGINI IPT2
  853. IPT2.ITYPEL = 22
  854. Nelem = 0
  855.  
  856. C Boucle sur les noeuds <<<<<<<<<<<<<<<
  857. DO 2012 INUM = 1, NBPTB
  858. IF (TBLOQ.MBLOQ(IENR, INUM).LT.2) goto 2012
  859. Nelem = Nelem + 1
  860. IPT2.NUM(2,Nelem)=INUM
  861. NBPTS1 = NBPTS1 + 1
  862. IPT2.NUM(1,Nelem)=NBPTS1
  863. 2012 CONTINUE
  864. C Fin de boucle sur les noeuds <<<<<<<<
  865.  
  866. C Matrice de blocage
  867. SEGINI XMATR2
  868. DO i=1,NELRIG
  869. XMATR2.RE(1,1,i)=0.D0
  870. XMATR2.RE(2,1,i)=1.D0
  871. XMATR2.RE(2,2,i)=0.D0
  872. XMATR2.RE(1,2,i)=1.D0
  873. ENDDO
  874. C Segment Descripteur
  875. SEGINI DES2
  876. DES2.LISINC(1)='LX'
  877. DES2.LISDUA(1)='FLX'
  878. DES2.LISINC(2)=MOTINC(IDDL)
  879. DES2.LISDUA(2)=MOTDUA(IDDL)
  880.  
  881. DES2. NOELEP(1)=1
  882. DES2. NOELEP(2)=2
  883. DES2. NOELED(1)=1
  884. DES2. NOELED(2)=2
  885.  
  886. C stockage de la rigidité construite dans LOCIRI
  887.  
  888. LOCIRI(1, NBRIG)=IPT2
  889. LOCIRI(3, NBRIG)=DES2
  890. LOCIRI(4, NBRIG)=XMATR2
  891. c desactivations
  892. SEGDES IPT2, DES2, XMATR2
  893.  
  894. 2022 CONTINUE
  895. IF (IDIM.EQ.2) IDDL=IDDL+1
  896. 2032 CONTINUE
  897. 2002 CONTINUE
  898. C write (*,*) '***************NBRIG*******************', NBRIG
  899. C Fin des boucles sur les niveau d'enrichissement et composantes <<<<<<<
  900.  
  901. NBPTS = NBPTS1
  902. SEGADJ MCOORD
  903. c
  904. C*********************************************************
  905. C SUPPRESSION ET DESACTIVATION DE SEGMENTS
  906. C*********************************************************
  907. C
  908. SEGSUP,WRK1,WRK2,MVELCH
  909. SEGSUP,MRACC
  910.  
  911. segdes,MELEME
  912. segdes,MINT1
  913. if(MINT2.ne.0) segdes,MINT2
  914. c
  915. c desactivation des maillages et segment imatri
  916. do ienr=1,(1+NBENR1)
  917. c IPT1 = LOCIRI(1,ienr)
  918. c if(IPT1.ne.0) segdes,IPT1
  919. XMATRI = LOCIRI(4,ienr)
  920. if(XMATRI.ne.0) segdes,XMATRI
  921. enddo
  922.  
  923. c SEGDES dans RIGI1 = IMODEL
  924. C
  925. RETURN
  926. END
  927.  
  928.  
  929.  
  930.  
  931.  
  932.  

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