Télécharger grad1x.eso

Retour à la liste

Numérotation des lignes :

grad1x
  1. C GRAD1X SOURCE OF166741 25/02/21 21:17:15 12166
  2. C
  3. subroutine GRAD1X (IMODEL,IVADEP,LRE,IVAGRA,NGRA,
  4. & MINT1,MINT2,IIPDPG,IOK)
  5. c
  6. c CALCUL DU GRADIENT DANS LE CAS D'ELEMENTS XFEM (MFR=63)
  7. c
  8. c largement inspiré de epsix.eso + bgrmas.eso
  9. C
  10. C
  11. C*********************************************************
  12. C PARTIE DECLARATIVE
  13. C*********************************************************
  14.  
  15. IMPLICIT INTEGER(I-N)
  16. IMPLICIT REAL*8(A-H,O-Z)
  17.  
  18. -INC PPARAM
  19. -INC CCOPTIO
  20. -INC SMCOORD
  21. -INC SMMODEL
  22. -INC SMCHAML
  23. -INC SMINTE
  24. -INC SMELEME
  25. -INC SMLREEL
  26.  
  27. POINTEUR MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE
  28.  
  29. -INC TMPTVAL
  30.  
  31. SEGMENT WRK1
  32. REAL*8 XE(3,NBBB)
  33. REAL*8 XDDL(LRE),GRADI(NGRA)
  34. ENDSEGMENT
  35. c
  36. SEGMENT WRK2
  37. REAL*8 SHPWRK(6,NBNO),BGR(NGRA,LRE)
  38. REAL*8 LV7WRK(NBENRMA2,2,6,NBBB)
  39. ENDSEGMENT
  40.  
  41. SEGMENT MRACC
  42. INTEGER TLREEL(NBENRMA2,NBI)
  43. ENDSEGMENT
  44.  
  45. PARAMETER (NDDLMAX=30,NBNIMAX=10)
  46. PARAMETER (NBENRMAX=5)
  47. DIMENSION MLRE(NBENRMAX+1)
  48.  
  49. C write (*,*) '############################'
  50. C write (*,*) '##### DEBUT DE GRAD1X #####'
  51. C write (*,*) '############################'
  52. C
  53. C*********************************************************
  54. c
  55. C RECUPERATION + ACTIVATIONS + VALEURS UTILES
  56. c
  57. C*********************************************************
  58. c
  59. IOK = 0
  60. C++++ Recup de la geometrie deja activée par grad1 ++++++
  61. MELEME= IMAMOD
  62. C nbre de noeuds par element, nbre d elements
  63. NBNN1 = NUM(/1)
  64. NBEL1 = NUM(/2)
  65.  
  66. C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++
  67. MELE = NEFMOD
  68. NGAU1 = MINT1.POIGAU(/1)
  69. c C sous decoupage et points de Gauss de l element geometrique de base
  70. c IF(MINT2.NE.0) SEGACT,MINT2
  71. c NGAU2 = MINT2.POIGAU(/1)
  72.  
  73. C++++ Recup des infos d enrichissement +++++++++++++++++++
  74. c recup du MCHEX1 d'enrichissement
  75. MCHAM1 = 0
  76. NOBMO1 = IVAMOD(/1)
  77. if (NOBMO1.ne.0) then
  78. do iobmo1=1,NOBMO1
  79. if ((TYMODE(iobmo1)).eq.'MCHAML') then
  80. MCHEX1 = IVAMOD(iobmo1)
  81. segact,MCHEX1
  82. if ((MCHEX1.TITCHE).eq.'ENRICHIS') then
  83. MCHAM1 = MCHEX1.ICHAML(1)
  84. segact,MCHAM1
  85. goto 1000
  86. endif
  87. endif
  88. enddo
  89. write(*,*) 'Le modele est vide (absence d enrichissement)'
  90. else
  91. write(*,*) 'Il n y a pas de MCHEML enrichissement dans le Modele'
  92. endif
  93.  
  94. 1000 continue
  95. c niveau d enrichissement(s) du modele (ddls std u exclus)
  96. c NBENR1= 0 si std, 1 si H seul, 2 si H+F1, 3 si H+F1+F2, etc...
  97. if (MCHAM1.ne.0) then
  98. NBENR1= MCHAM1.IELVAL(/1)
  99. else
  100. NBENR1 = 0
  101. endif
  102. C write(*,*) 'niveau d enrichissement(s) du modele',NBENR1
  103. C
  104. C*********************************************************
  105. C INITIALISATIONS...
  106. C*********************************************************
  107. c
  108. c preparation des tables avec:
  109.  
  110. if (NBENR1.ne.0) then
  111. do ienr=1,NBENR1
  112. c -le nombre d'inconnues de chaque sous-zone
  113. c determinee depuis le nombre de fonction de forme
  114. c ienr= 1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10)
  115. nbniJ = 2 + ((ienr-1)*4)
  116. MLRE(1+ienr) = IDIM*NBNN1*nbniJ
  117. enddo
  118. endif
  119. C Tables + longues car 1er indice correspond au fontion de forme std
  120. MLRE(1) = IDIM*NBNN1*1
  121.  
  122. if (NBENR1.lt.(NBENRMAX+1)) then
  123. do ienr=(NBENR1+1),(NBENRMAX)
  124. MLRE(1+ienr) = 0
  125. enddo
  126. endif
  127. c
  128. c ...DU SEGMENT WRK1
  129. NBENRMA2 = NBENRMAX
  130. NBBB = NBNN1
  131. SEGINI,WRK1
  132.  
  133. c ...DU SEGMENT WRK2
  134. c NBNO = NBNI
  135. NBNO = LRE/IDIM
  136. SEGINI,WRK2
  137. C
  138. c ...DU SEGMENT MRACC
  139. NBENRMA2 = NBENRMAX
  140. NBI = NBNN1
  141. segini,MRACC
  142. C
  143. C du nombre d erreur sur le noms de composantes
  144. NBERR1=0
  145.  
  146.  
  147.  
  148. C*********************************************************
  149. C
  150. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BOUCLE SUR LES ELEMENTS
  151. C
  152. DO 2000 J=1,NBEL1
  153. C write(*,*) '========= EF',J,' /',NBEL1,'========='
  154. C
  155. C
  156. C*********************************************************
  157. C POUR CHAQUE ELEMENT, ...
  158. C*********************************************************
  159. C
  160. C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT
  161. CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE)
  162. C
  163. c
  164. C++++ NBENRJ = niveau d enrichissement de l element ++++
  165. C =0 si EF std =1 si U+H =2 si U+H+F1 =3 si U+H+F1+F2
  166. NBENRJ=0
  167. if(NBENR1.ne.0) then
  168. do IENR=1,NBENR1
  169. MELVA1 = MCHAM1.IELVAL(IENR)
  170. segact,MELVA1
  171. N2PTEL=MELVA1.IELCHE(/1)
  172. N2EL =MELVA1.IELCHE(/2)
  173. do I=1,NBNN1
  174. IF (N2PTEL.EQ.1 .AND. N2EL.EQ.1)THEN
  175. C Cas du champ constant sur tout le MAILLAGE support
  176. mlree1 = MELVA1.IELCHE(1,1)
  177. ELSEIF(N2PTEL.EQ.1)THEN
  178. C Cas du champ constant par element
  179. mlree1 = MELVA1.IELCHE(1,J)
  180. ELSE
  181. C Cas du champ variable sur les noeuds et les elements
  182. mlree1 = MELVA1.IELCHE(I,J)
  183. ENDIF
  184.  
  185. C on en profite pour remplir MRACC table de raccourcis pour cet element
  186. TLREEL(IENR,I) = mlree1
  187. if (mlree1.ne.0) then
  188. NBENRJ = max(NBENRJ,IENR)
  189. C et on active la listreel
  190. segact,mlree1
  191. endif
  192. enddo
  193. enddo
  194. endif
  195. if (NBENRMA2.gt.NBENR1) then
  196. do IENR2=(NBENR1+1),NBENRMA2
  197. do I=1,NBNN1
  198. TLREEL(IENR2,I) = 0
  199. enddo
  200. enddo
  201. endif
  202. C
  203. c
  204. c++++ quelques indices pour dimensionner au plus juste
  205. c nbre total de ddl de l'élément considéré
  206. c NLIGRD = MLRE(1+NBENRJ)
  207. c NLIGRD = MLRE(1+NBENRJ)
  208. NDDL = MLRE(1+NBENRJ)
  209. c nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim
  210. NBNI = (MLRE(1+NBENRJ)) / IDIM
  211. * nbre d inconnue / noeud
  212. NCOMP = NDDL/NBNN1
  213. C write(6,*) 'EF',J,' NBENRJ=',NBENRJ,
  214. C &'donc',NDDL,' ddls et ',NBNI,' fonctions de forme'
  215. c
  216. c
  217. C
  218. C++++ ON CALCULE XDDL ++++
  219. MPTVAL = IVADEP
  220. NCOSOU = IVAL(/1)
  221. C
  222. INITOT = 0
  223. C-->> BOUCLE SUR LES niveaux d'ENRICHISSEMENTS (0:U 1:A 2:BCDE1 3:BCDE2)
  224. DO IENR=0,NBENRJ
  225. *nbre de fonction(s) de ce niveau d'enrichissement (=1 si std ou H, =4 pour F1,2,...)
  226. if (IENR .le. 1) then
  227. NBNIENR = 1
  228. else
  229. NBNIENR = 4
  230. endif
  231. C---->> BOUCLE SUR LES fonctions de forme de ce type d'enrichissement
  232. DO INI=1,NBNIENR
  233. INITOT = INITOT + 1
  234. C------>> BOUCLE SUR LA DIMENSION
  235. DO 2220 KDIM=1,IDIM
  236. ICOMP = (INITOT-1)*IDIM + KDIM
  237.  
  238. c --cas ou on n'a pas trouvé assez de composantes--
  239. if(ICOMP.GT.NCOSOU) GOTO 2221
  240.  
  241. MELVAL = IVAL(ICOMP)
  242. c --cas ou on a pas trouvé cette composante dans cette zone du
  243. c chpoint solution devenu mchaml --
  244. if(MELVAL.eq.0)then
  245. c Avait on besoin de cette composante?
  246. c oui, si c'est une composante obligatoire
  247. if(IENR.eq.0) goto 2991
  248. c oui, si l'un des noeuds est enrichi
  249. do I=1,NBNN1
  250. if(TLREEL(IENR,I).ne.0) goto 2991
  251. enddo
  252. c non, si c'est facultatif et qu'on n'est pas enrichi -> on saute
  253. goto 2220
  254. c ->AVERTISSEMENT puis on saute
  255. 2991 NBERR1=NBERR1+1
  256. if(IIMPI.lt.1) goto 2220
  257. c write(IOIMP,*) 'PB OPERATEUR GRAD :'
  258. write(IOIMP,991) ICOMP,IENR,INI,KDIM
  259. 991 format(2X,'ABSENCE DANS LE CHPOINT DEPLACEMENT DE LA ',I3,
  260. $ 'ieme COMPOSANTE (enrichissement',I3,
  261. $ ', fonction',I3,', direction ',I3,
  262. $ ') NECESSAIRE POUR L UN DES NOEUDS SUIVANTS :')
  263. write(IOIMP,*)' noeuds :',(NUM(iou,J),iou=1,NBNN1)
  264. goto 2220
  265. endif
  266.  
  267. C---------->> BOUCLE SUR LES NOEUDS
  268. N1PTEL = VELCHE(/1)
  269. N1EL = VELCHE(/2)
  270. DO I=1,NBNN1
  271. IQ = ((INITOT-1)*NBNN1*IDIM) + ((I-1)*IDIM) + KDIM
  272. IF (N1PTEL.EQ.1 .AND. N1EL.EQ.1)THEN
  273. C Cas du champ constant sur tout le MAILLAGE support
  274. XDDL(IQ) = VELCHE(1,1)
  275. ELSEIF(N1PTEL.EQ.1) THEN
  276. C Cas du champ constant par element
  277. XDDL(IQ) = VELCHE(1,J)
  278. ELSE
  279. C Cas du champ variable sur les noeuds et les elements
  280. XDDL(IQ) = VELCHE(I,J)
  281. ENDIF
  282. ENDDO
  283.  
  284. 2220 CONTINUE
  285. ENDDO
  286. ENDDO
  287.  
  288. c --cas normal (toutes les composantes souhaitees etaient presentes)--
  289. GOTO 2223
  290.  
  291. c --cas ou on n'a pas trouvé assez de composantes--
  292. 2221 CONTINUE
  293. if (IIMPI.ge.1) then
  294. WRITE(IOIMP,2222) J,NCOSOU,NCOMP
  295. 2222 FORMAT('OPERATEUR GRAD : ELEMENT ',I6,
  296. & ' LE CHAMP DE DEPLACEMENT CONTIENT ',I3,' COMPOSANTES',
  297. & ' PAR NOEUD AU LIEU DE ',I3,' ATTENDUES')
  298. endif
  299. NDDL=NCOSOU*NBNN1
  300. NBENRJ=IENR
  301.  
  302. 2223 CONTINUE
  303. c --cas ou on a une ou des erreurs--
  304. IF (NBERR1.gt.0.and.J.eq.NBEL1) THEN
  305. write(IOIMP,*) 'OPERATEUR GRAD : ABSENCE DANS LE CHPOINT ',
  306. & 'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE'
  307. ENDIF
  308.  
  309. c
  310. c
  311. c rem: il serait probablement interessant au niveau du tems cpu
  312. c d'utiliser moins de pts de Gauss lorsque l element est élastique
  313. c On pourrait par ex. utiliser MINT2 = infell(12) qui contient
  314. c le segment d'integration de l'EF std (QUA4 par ex.)
  315. * if((NBENRJ.eq.0).and.(MINT2.ne.0)) then
  316. * MINTE = MINT2
  317. * NBPGAU= NGAU2
  318. * else
  319. MINTE = MINT1
  320. NBPGAU= NGAU1
  321. * endif
  322. c
  323. C
  324. C*********************************************************
  325. C
  326. C>>>>>>>>>> BOUCLE SUR LES POINTS DE GAUSS
  327. C
  328. DO 2500 KGAU=1,NBPGAU
  329. C
  330. C*********************************************************
  331. C Initialisation à 0
  332. C*********************************************************
  333.  
  334. c ZERO ne serait-il pas facultatif?
  335. CALL ZERO(SHPWRK,6,NBNO)
  336. XGAU = 0.D0
  337. YGAU = 0.D0
  338. ZGAU = 0.D0
  339. C
  340. i6zz = 3
  341. IF(IDIM.EQ.3) i6zz = 4
  342. C
  343. do ienr7=1,NBENRJ
  344. do inod7=1,NBNN1
  345. do i6=1,i6zz
  346. LV7WRK(ienr7,1,i6,inod7) = 0.D0
  347. LV7WRK(ienr7,2,i6,inod7) = 0.D0
  348. enddo
  349. enddo
  350. enddo
  351.  
  352.  
  353. c*********************************************************
  354. c Calcul des fonction de forme std dans repere local
  355. c*********************************************************
  356.  
  357. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  358. c (et donc sur les Ni std)
  359. DO 2510 I=1,NBNN1
  360.  
  361. C++++ Calcul des Ni std
  362. c (rappel: 1:Ni 2:Ni,qsi 3:Ni,eta avec i=1,4)
  363. SHPWRK(1,I) = SHPTOT(1,I,KGAU)
  364. SHPWRK(2,I) = SHPTOT(2,I,KGAU)
  365. SHPWRK(3,I) = SHPTOT(3,I,KGAU)
  366. IF(IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)
  367.  
  368. C++++ Calcul des coordonnees dans le repere global
  369. XGAU = XGAU + (SHPWRK(1,I)*XE(1,I))
  370. YGAU = YGAU + (SHPWRK(1,I)*XE(2,I))
  371. IF(IDIM.EQ.3) ZGAU = ZGAU + (SHPWRK(1,I)*XE(3,I))
  372.  
  373. 2510 CONTINUE
  374. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  375. c if(J.eq.77) write(6,*) 'xgau_',KGAU,' =',XGAU,YGAU,ZGAU
  376.  
  377.  
  378.  
  379. c*********************************************************
  380. c Passage des fonctions de forme std dans repere global
  381. c*********************************************************
  382.  
  383. C++++ CALCUL DES Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  384. CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  385. c if(J.eq.1.and.KGAU.eq.1)
  386. c & write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1)
  387.  
  388. c*********************************************************
  389. c Si on est pas du tout enrichi on peut sauter une grosse partie
  390. c*********************************************************
  391. if(NBENRJ.eq.0) goto 2999
  392.  
  393. c*********************************************************
  394. c Calcul des level set + leurs derivees dans repere global
  395. c*********************************************************
  396.  
  397. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  398. do 2520 ienr=1,NBENRJ
  399.  
  400. c MELVA1=MCHAM1.IELVAL(IENR)
  401. c segact,MELVA1
  402.  
  403. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  404. DO 2521 I=1,NBNN1
  405.  
  406. C++++ Le I eme noeud est-il ienr-enrichi?
  407. mlree1= TLREEL(IENR,I)
  408. c if(J.eq.77) write(6,*) 'ienr,I',ienr,I,' mlree1=',mlree1
  409. if(mlree1.eq.0) goto 2521
  410.  
  411.  
  412. c Calcul du repere local de fissure(=PSI,PHI)
  413. c (rappel: 1,1:psi 1,2:phi 2,1 psi,x ...etc...)
  414. do 2522 inode=1,NBNN1
  415. c pour le H-enrichissement, on n a pas gardé PSI (inutile)
  416. if (ienr.ne.1) then
  417. c valeur de PSI au inode^ieme noeud
  418. xpsi1 = mlree1.PROG(inode)
  419. c qu on multiplie par la valeur de Ni^std au pt de G considéré
  420. LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
  421. & + (SHPWRK(1,inode)*xpsi1)
  422. LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
  423. & + (SHPWRK(2,inode)*xpsi1)
  424. LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
  425. & + (SHPWRK(3,inode)*xpsi1)
  426. IF(IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  427. & + (SHPWRK(4,inode)*xpsi1)
  428. c valeur de PHI au inode^ieme noeud
  429. xphi1 = mlree1.PROG(NBNN1+inode)
  430. else
  431. xphi1 = mlree1.PROG(inode)
  432. endif
  433. LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
  434. & + (SHPWRK(1,inode)*xphi1)
  435. LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
  436. & + (SHPWRK(2,inode)*xphi1)
  437. LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
  438. & + (SHPWRK(3,inode)*xphi1)
  439. IF(IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  440. & + (SHPWRK(4,inode)*xphi1)
  441. 2522 continue
  442.  
  443. 2521 continue
  444. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  445.  
  446.  
  447. 2520 CONTINUE
  448. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  449.  
  450. c on a construit
  451. C LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  452.  
  453.  
  454. c*********************************************************
  455. c Ajout des fonctions de forme d enrichissement
  456. c + leurs derivees dans repere global
  457. c*********************************************************
  458. CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET)
  459.  
  460.  
  461.  
  462. c retour a la partie commune aux EF enrichi et non enrichi
  463. 2999 continue
  464.  
  465. C*********************************************************
  466. C CALCUL DE BGR
  467. C*********************************************************
  468. c ZERO ne serait-il pas facultatif?
  469. call ZERO(BGR,9,NDDL)
  470. KB=1
  471. C boucle sur tous les Ni
  472. DO 3001 II=1,NBNI
  473.  
  474. c partie commune 2D 3D
  475. BGR(1,KB) = SHPWRK(2,II)
  476. BGR(2,KB) = SHPWRK(3,II)
  477. BGR(4,KB+1) = SHPWRK(2,II)
  478. BGR(5,KB+1) = SHPWRK(3,II)
  479.  
  480. c cas 3D
  481. IF (IDIM.EQ.3) THEN
  482. BGR(3,KB) = SHPWRK(4,II)
  483. BGR(6,KB+1) = SHPWRK(4,II)
  484. BGR(7,KB+2) = SHPWRK(2,II)
  485. BGR(8,KB+2) = SHPWRK(3,II)
  486. BGR(9,KB+2) = SHPWRK(4,II)
  487. ENDIF
  488.  
  489. c cas 2D PLAN GENE
  490. IF (IFOUR.EQ.-3) THEN
  491. IREF=(IIPDPG-1)*(IDIM+1)
  492. BGR(9,KB) = 1.D0
  493. BGR(9,KB+1) = XCOOR(IREF+1)-XGAU
  494. BGR(9,KB+2) = YGAU-XCOOR(IREF+2)
  495. ENDIF
  496.  
  497. KB = KB + IDIM
  498.  
  499. 3001 CONTINUE
  500. C
  501. c
  502.  
  503. C*********************************************************
  504. C CALCUL DE BGR.q
  505. C*********************************************************
  506.  
  507. C APPEL A LA PROCEDURE DE CALCUL
  508. CALL BGRDEP(BGR,NGRA,XDDL,NDDL,GRADI)
  509. c
  510.  
  511. C*********************************************************
  512. C ECRITURE DES DIFFERENTES COMPOSANTES DU GRADIENT
  513. C*********************************************************
  514. MPTVAL=IVAGRA
  515. DO i=1,NGRA
  516. MELVAL=IVAL(i)
  517. IGMN=MIN(KGAU,VELCHE(/1))
  518. IBMN=MIN(J,VELCHE(/2))
  519. VELCHE(IGMN,IBMN)=GRADI(i)
  520. ENDDO
  521. c
  522. 2500 CONTINUE
  523. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  524. C
  525. 2000 CONTINUE
  526. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  527.  
  528. C*********************************************************
  529. C SUPPRESSION ET DESACTIVATION DE SEGMENTS
  530. C*********************************************************
  531. SEGSUP,WRK1,WRK2
  532. SEGSUP,MRACC
  533.  
  534. c on met mint2 a zero pour eviter pb dans grad1
  535. MINT2 = 0
  536. c
  537. c pas de retour en erreur pour linstant
  538. IOK=1
  539. c write(6,*) 'iok=',IOK
  540. C
  541. RETURN
  542. END
  543.  
  544.  
  545.  

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