Télécharger piocax.eso

Retour à la liste

Numérotation des lignes :

  1. C PIOCAX SOURCE PV 18/11/20 21:15:04 1001
  2. SUBROUTINE PIOCAX(NBNN,IDIM,TAB1,NCOELE,NBPTEL,IPMINT,XE1,XE2,
  3. 1 TABA,MRACC,SH1,TAB,MWRK6,LHOOK,
  4. 2 IFOU,KCAS,KERRE)
  5. C=======================================================================
  6. C
  7. C TRANSFORME LES CONTRAINTES DE PIOLA KIRCHHOFF EN CONTRAINTES DE
  8. C CAUCHY
  9. C ENTREE
  10. C -------
  11. C NBNN = NOMBRE DE POINTS PAR ELEMENTS
  12. C IDIM = DIMENSION DE L ESPACE SUPPORT
  13. C
  14. C TAB1(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE PIOLA KIRCHHOFF
  15. C
  16. C NCOELE = NOMBRE DE COMPOSTS TABLEAU DES CONTRAINTES
  17. C
  18. C NBPTEL = NOMBRE DE POINTS DE GAUSS
  19. C IPMINT = POINTEUR DES FONCTIONS DE FORME
  20. C TABA = pointeur tableau avec ddl de saut
  21. c MRACC = pointeur tableau de raccourci pour les
  22. C enrichissements elementaires
  23. C
  24. C KCAS = 1 SI CONTRAINTES, 2 SI DEFORMATIONS
  25. C
  26. C TABLEAUX DE TRAVAIL
  27. C--------------------
  28. C XE1(3,NBNN) = COORDONNEES CORRESPONDANT A LA CONFIGURATION DEPART
  29. C
  30. C XE2(3,NBNN) = COORDONNEES CORRESPONDANT A LA CONFIGURATION ACTUEL
  31. C
  32. C SH1(6,NBNN) = FONCTIONS DE FORME EN UN POINT DE GAUSS
  33. C
  34. C SORTIES
  35. C---------
  36. C TAB(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE CAUCHY
  37. C
  38. C
  39. C AOUT 85
  40. C MODIF PEGON FEV 90 CAS BIDIM
  41. C PASSAGE AUX NOUVEAUX CHAMELEMS PAR P.DOWLATYARI 12/4/91
  42. C
  43. C=======================================================================
  44. IMPLICIT INTEGER(I-N)
  45. IMPLICIT REAL*8(A-H,O-Z)
  46. C
  47. -INC CCREEL
  48. -INC SMLREEL
  49. -INC SMINTE
  50. *
  51. SEGMENT MWRK6
  52. INTEGER ITRES1(NBPTEL)
  53. REAL*8 PRODDI(NBPTEL,LHOO2),PRODDO(NBPTEL,LHOO2)
  54. REAL*8 DDHOOK(LHOOK,LHOOK),DDHOMU(LHOOK,LHOOK)
  55. REAL*8 VEC(LHOOK),VEC2(LHOOK)
  56. ENDSEGMENT
  57. *
  58.  
  59. DIMENSION TAB1(NBPTEL,*)
  60. DIMENSION TAB(NBPTEL,*)
  61. DIMENSION XE1(3,*),XE2(3,*),SH1(6,*)
  62. *as xfem 2010_01_13
  63. PARAMETER (NBNNMAX=20)
  64. SEGMENT MRACC
  65. INTEGER TLREEL(NBNN)
  66. ENDSEGMENT
  67. * ddl de saut enrichissement
  68. SEGMENT TABA
  69. REAL*8 TABA1(IDIM,NBNN),TABA2(IDIM,NBNN)
  70. ENDSEGMENT
  71.  
  72. * phi et H aux nnoeuds ; phi et H au point de Gauss courant
  73. DIMENSION xphi1(NBNNMAX),xh1(NBNNMAX),TABPHG(2),tabdh(nbnnMAX)
  74. *fin as xfem 2010_01_13
  75. C
  76. C TABLEAUX DE TRAVAIL DIMENSIONNES ICI
  77. C
  78. DIMENSION XJAC(3,3),FAC(6)
  79. C
  80. C TABLEAUX INDIQUANT LA CORRESPONDANCE ENTRE INDICES I,J ET NUMERO
  81. C DE LA COMPOSANTE DE CONTRAINTES OU DE DEFORMATIONS
  82. C
  83. DIMENSION IN(6),JN(6),ITAB(3,3)
  84. C
  85. DATA FAC/1.D0,1.D0,1.D0,0.5D0,0.5D0,0.5D0/
  86. DATA IN/1,2,3,1,1,2/
  87. DATA JN/1,2,3,2,3,3/
  88. C
  89. DATA ITAB(1,1),ITAB(1,2),ITAB(1,3)/1,4,5/
  90. DATA ITAB(2,1),ITAB(2,2),ITAB(2,3)/4,2,6/
  91. DATA ITAB(3,1),ITAB(3,2),ITAB(3,3)/5,6,3/
  92. C
  93. KERRE=0
  94.  
  95. MINTE = IPMINT
  96. NBSH = SHPTOT(/2)
  97.  
  98. IDIM1=IDIM+1
  99. C
  100. C MISE A ZERO DES CONTRAINTES OU DES DEFORMATIONS
  101. C
  102. DO 50 IB=1,NCOELE
  103. DO 50 IA=1,NBPTEL
  104. TAB(IA,IB)=0.D0
  105. 50 CONTINUE
  106. C
  107. C BOUCLE SUR LES POINTS DE GAUSS
  108. C
  109. DO 130 IC=1,NBPTEL
  110.  
  111. tabphg(2)=0.D0
  112. *as xfem 2010_01_13
  113. * Initialisation de SH1(Ni, Ni,x, Ni,y)
  114. do i3=1,nbnn
  115. xh1(i3)=0.D0
  116. do i4=1,IDIM1
  117. SH1(i4,i3)=SHPTOT(i4,i3,IC)
  118. enddo
  119. enddo
  120. * Calcul de H et phi aux noeuds et au point de Gauss IC
  121. do 131 i1=1,nbnn
  122. mlree1=tlreel(i1)
  123. if(mlree1.eq.0) goto 131
  124. tabphg(1)=0.D0
  125. do i2=1,nbnn
  126. xphi1(i2)=mlree1.PROG(i2)
  127. if (abs(xphi1(i2)).lt.1.d-7) then
  128. xh1(i2)=0.D0
  129. else
  130. xh1(i2)=sign(1.d0,xphi1(i2))
  131. endif
  132. tabphg(1)=tabphg(1)+SH1(1,i2)*xphi1(i2)
  133. enddo
  134. if (abs(tabphg(1)).lt.1.d-7) then
  135. tabphg(2)=0.D0
  136. else
  137. tabphg(2)=sign(1.d0,tabphg(1))
  138. endif
  139.  
  140. 131 continue
  141. * Calcul des H(x)-H(xi) :
  142. do i3=1,nbnn
  143. tabdh(i3)=tabphg(2)-xh1(i3)
  144. enddo
  145. * Calcul de SH1 :
  146. call jacobix(XE1,TABA1,TABDH,SH1,IDIM,NBNN,DJac)
  147. C
  148. C CALCUL DE LA MATRICE F
  149. C
  150. CALL ZERO(XJAC,3,3)
  151. DO 140 ID=1,NBNN
  152. DO 140 IE=1,IDIM
  153. * r_z = XE2(IE,ID)
  154. r_z = XE2(IE,ID)+(tabdh(ID)*TABA2(IE,ID))
  155. DO 140 IF=1,IDIM
  156. XJAC(IE,IF)=XJAC(IE,IF) + SH1(IF+1,ID)*r_z
  157. 140 CONTINUE
  158. *fin as xfem 2010_01_13
  159. IF(IDIM.EQ.2) THEN
  160. XJAC(3,3)=1.D0
  161. IF(IFOU.EQ.0) THEN
  162. C
  163. CCCCCCCCCCCCC CAS AXISYMETRIQUE
  164. C
  165. R1=0.
  166. R2=0.
  167. DO 150 ID=1,NBNN
  168. R1=R1+SH1(1,ID)*XE1(1,ID)
  169. R2=R2+SH1(1,ID)*XE2(1,ID)
  170. 150 CONTINUE
  171. if (r1.lt.-xpetit/xzprec.or.r1.gt.xpetit/xzprec)then
  172. XJAC(3,3)=R2/R1
  173. else
  174. xjac(3,3)=xgrand*xzprec
  175. endif
  176. ENDIF
  177. ENDIF
  178. C
  179. C
  180. GO TO (500,600,700),KCAS
  181. C
  182. C KCAS=1 CAS DES CONTRAINTES
  183. C ----------------------------
  184. C
  185. 500 CONTINUE
  186. C
  187. C
  188. CCCCCCCCCCCC CALCUL DE DETERMINANT DE F
  189. C
  190. IF(IDIM.EQ.2) THEN
  191. DETF=XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1)
  192. DETF = DETF * XJAC (3,3)
  193. ENDIF
  194. IF(IDIM.EQ.3) THEN
  195. DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  196. DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  197. DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  198. ENDIF
  199. if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  200. DETF=1./(DETF)
  201. else
  202. DETF=XGRAND*xzprec
  203. endif
  204. C
  205. C CALCUL DES CONTRAINTES DE CAUCHY
  206. C
  207. DO 160 ID=1,NCOELE
  208. IND=IN(ID)
  209. JND=JN(ID)
  210. DO 170 IE=1,IDIM
  211. DO 170 IF=1,IDIM
  212. ICO=ITAB(IE,IF)
  213. TAB(IC,ID)=TAB1(IC,ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF
  214. 1 +TAB(IC,ID)
  215. 170 CONTINUE
  216. 160 CONTINUE
  217. C
  218. C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  219. C
  220. IF(IDIM.EQ.2) THEN
  221. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  222. ENDIF
  223. GO TO 130
  224. C
  225. C KCAS=2 CAS DES DEFORMATIONS
  226. C -----------------------------
  227. C
  228. 600 CONTINUE
  229. C
  230. C
  231. CCCCCCCCCCCC CALCUL DE L'INVERSE DE F
  232. C
  233. CALL INVMA1(XJAC,3,3,KERRE)
  234. IF(KERRE.NE.0) THEN
  235. WRITE(6,77881) ((XJAC(MI,MJ),MJ=1,3),MI=1,3)
  236. 77881 FORMAT(2X,' MATRICE SINGULIERE' /(3(1X,1PE12.5)/))
  237. RETURN
  238. ENDIF
  239. C
  240. C CALCUL DES DEFORMATIONS
  241. C
  242. DO 260 ID=1,NCOELE
  243. IND=IN(ID)
  244. JND=JN(ID)
  245. DO 270 IE=1,IDIM
  246. DO 270 IF=1,IDIM
  247. ICO=ITAB(IE,IF)
  248. TAB(IC,ID)=TAB(IC,ID) +
  249. 1 FAC(ICO)*TAB1(IC,ICO)*XJAC(IE,IND)*XJAC(IF,JND)/FAC(ID)
  250. 270 CONTINUE
  251. 260 CONTINUE
  252. C
  253. C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  254. C
  255. IF(IDIM.EQ.2) THEN
  256. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)
  257. ENDIF
  258. C
  259. GO TO 130
  260. C
  261. C KCAS=3 CAS DE LA MATRICE DE HOOKE
  262. C -----------------------------------
  263. C
  264. 700 CONTINUE
  265. C
  266. CCCCCCCCCCCC CALCUL DE L'INVERSE DU DETERMINANT DE F
  267. C
  268. IF (IDIM.EQ.3) THEN
  269. DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  270. DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  271. DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  272. ELSE IF (IDIM.EQ.2) THEN
  273. DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3)
  274. ELSE IF (IDIM.EQ.1) THEN
  275. DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3)
  276. ENDIF
  277. if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  278. DETF=1./(DETF)
  279. else
  280. DETF=XGRAND*xzprec
  281. endif
  282. C
  283. IJ=1
  284. DO 710 JJ=1,LHOOK
  285. DO 710 II=1,LHOOK
  286. DDHOOK(II,JJ)=PRODDI(IC,IJ)
  287. IJ=IJ+1
  288. 710 CONTINUE
  289. *
  290. CALL ZERO(DDHOMU,LHOOK,LHOOK)
  291. KEY=2
  292. C
  293. DO 720 LC=1,LHOOK
  294.  
  295. CALL ZERO (VEC,LHOOK,1)
  296. DO 760 ID=1,LHOOK
  297. IND=IN(ID)
  298. JND=JN(ID)
  299. DO 770 IE=1,IDIM
  300. DO 770 IF=1,IDIM
  301. ICO=ITAB(IE,IF)
  302. IF(ICO.EQ.LC) THEN
  303. VEC(ID)=VEC(ID)+
  304. & FAC(ICO)*XJAC(IE,IND)*XJAC(IF,JND)/FAC(ID)
  305. ENDIF
  306. 770 CONTINUE
  307. 760 CONTINUE
  308. C
  309. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  310. C
  311. IF (IDIM.EQ.2) THEN
  312. IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3)
  313. *
  314. ELSE IF (IDIM.EQ.1) THEN
  315. IF(LC.EQ.2) VEC(2)=XJAC(2,2)*XJAC(2,2)
  316. IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3)
  317. ENDIF
  318. C
  319. CALL MATVE1(DDHOOK,VEC,LHOOK,LHOOK,VEC2,KEY)
  320. C
  321. DO 761 ID=1,LHOOK
  322. IND=IN(ID)
  323. JND=JN(ID)
  324. DO 771 IE=1,IDIM
  325. DO 771 IF=1,IDIM
  326. ICO=ITAB(IE,IF)
  327. DDHOMU(ID,LC)=VEC2(ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF
  328. 1 +DDHOMU(ID,LC)
  329. 771 CONTINUE
  330. 761 CONTINUE
  331. C
  332. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  333. C
  334. IF (IDIM.EQ.2) THEN
  335. DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  336. ELSE IF (IDIM.EQ.1) THEN
  337. DDHOMU(2,LC)=VEC2(2)*XJAC(2,2)*XJAC(2,2)*DETF
  338. DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  339. ENDIF
  340.  
  341. 720 CONTINUE
  342. C
  343. IJ=1
  344. DO 780 JJ=1,LHOOK
  345. DO 780 II=1,LHOOK
  346. PRODDO(IC,IJ)=DDHOMU(II,JJ)
  347. IJ=IJ+1
  348. 780 CONTINUE
  349. *
  350. *
  351. 130 CONTINUE
  352. RETURN
  353. END
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  

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