Télécharger piocaf.eso

Retour à la liste

Numérotation des lignes :

  1. C PIOCAF SOURCE PV 20/09/14 21:15:04 10713
  2. SUBROUTINE PIOCAF(NBNN,nbsh,IDIM,TAB1,NCOELE,NBPTEL,SHP,XE1,
  3. 1 XE2,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 NBSH = NOMBRE DE fonctions de forme
  13. C IDIM = DIMENSION DE L ESPACE SUPPORT
  14. C
  15. C TAB1(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE PIOLA KIRCHHOFF
  16. C
  17. C NCOELE = NOMBRE DE COMPOSTS TABLEAU DES CONTRAINTES
  18. C
  19. C NBPTEL = NOMBRE DE POINTS DE GAUSS
  20. C SHP(6,NBNN,NBPTEL)= FONCTIONS DE FORME
  21. C
  22. C KCAS = 1 SI CONTRAINTES, 2 SI DEFORMATIONS
  23. C
  24. C TABLEAUX DE TRAVAIL
  25. C--------------------
  26. C XE1(3,NBNN) = COORDONNEES CORRESPONDANT A LA CONFIGURATION DEPART
  27. C
  28. C XE2(3,NBNN) = COORDONNEES CORRESPONDANT A LA CONFIGURATION ACTUEL
  29. C
  30. C SH1(6,NBNN) = FONCTIONS DE FORME EN UN POINT DE GAUSS
  31. C
  32. C SORTIES
  33. C---------
  34. C TAB(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE CAUCHY
  35. C
  36. C
  37. C AOUT 85
  38. C MODIF PEGON FEV 90 CAS BIDIM
  39. C PASSAGE AUX NOUVEAUX CHAMELEMS PAR P.DOWLATYARI 12/4/91
  40. C
  41. C=======================================================================
  42. IMPLICIT INTEGER(I-N)
  43. IMPLICIT REAL*8(A-H,O-Z)
  44. -INC CCREEL
  45. *
  46. SEGMENT MWRK6
  47. INTEGER ITRES1(NBPTEL)
  48. REAL*8 PRODDI(NBPTEL,LHOO2),PRODDO(NBPTEL,LHOO2)
  49. REAL*8 DDHOOK(LHOOK,LHOOK),DDHOMU(LHOOK,LHOOK)
  50. REAL*8 VEC(LHOOK),VEC2(LHOOK)
  51. ENDSEGMENT
  52. *
  53. C
  54. dimension r_z(3)
  55. DIMENSION TAB1(NBPTEL,*)
  56. DIMENSION TAB(NBPTEL,*)
  57. *as xfem 2010_01_13
  58. DIMENSION SHP(6,NBSH,*)
  59. * DIMENSION SHP(6,NBNN,*)
  60. *fin as xfem 2010_01_13
  61. DIMENSION XE1(3,*),XE2(3,*),SH1(6,*)
  62. C
  63. C TABLEAUX DE TRAVAIL DIMENSIONNES ICI
  64. C
  65. DIMENSION XJAC(3,3),FAC(6),FACI(6)
  66. C
  67. C TABLEAUX INDIQUANT LA CORRESPONDANCE ENTRE INDICES I,J ET NUMERO
  68. C DE LA COMPOSANTE DE CONTRAINTES OU DE DEFORMATIONS
  69. C
  70. DIMENSION IN(6),JN(6),ITAB(3,3)
  71. C
  72. DATA FAC/1.D0,1.D0,1.D0,0.5D0,0.5D0,0.5D0/
  73. DATA FACI/1.D0,1.D0,1.D0,2.D0,2.D0,2.D0/
  74. DATA IN/1,2,3,1,1,2/
  75. DATA JN/1,2,3,2,3,3/
  76. C
  77. DATA ITAB(1,1),ITAB(2,1),ITAB(3,1)/1,4,5/
  78. DATA ITAB(1,2),ITAB(2,2),ITAB(3,2)/4,2,6/
  79. DATA ITAB(1,3),ITAB(2,3),ITAB(3,3)/5,6,3/
  80. C
  81. KERRE=0
  82. C
  83. C MISE A ZERO DES CONTRAINTES OU DES DEFORMATIONS
  84. C
  85. ** DO 50 IA=1,NBPTEL
  86. ** DO 50 IB=1,NCOELE
  87. ** TAB(IA,IB)=0.D0
  88. **50 CONTINUE
  89. CALL ZERO(TAB,NBPTEL,NCOELE)
  90. C
  91. C BOUCLE SUR LES POINTS DE GAUSS
  92. C
  93. DO 130 IC=1,NBPTEL
  94. C
  95. if (nbsh.ne.nbnn) then
  96. CALL HPRIMEX(XE1,NBNN,IDIM,NBSH,SHP,IC,SH1,DJAC)
  97. else
  98. CALL HPRIME(XE1,NBNN,IDIM,SHP,IC,SH1,DJAC)
  99. endif
  100. C
  101. C CALCUL DE LA MATRICE F (=XJAC)
  102. C
  103. CALL ZERO(XJAC,3,3)
  104. DO IF=1,IDIM
  105. JF = IF + 1
  106. r_z1 = 0.D0
  107. r_z2 = 0.D0
  108. r_z3 = 0.D0
  109. do id=1,nbnn
  110. r_z1 = r_z1 + SH1(JF,ID)*XE2(1,ID)
  111. if (idim.ge.2) r_z2 = r_z2 + SH1(JF,ID)*XE2(2,ID)
  112. if (idim.ge.3) r_z3 = r_z3 + SH1(JF,ID)*XE2(3,ID)
  113. enddo
  114. if (idim.ge.1) XJAC(1,IF) = r_z1
  115. if (idim.ge.2) XJAC(2,IF) = r_z2
  116. if (idim.ge.3) XJAC(3,IF) = r_z3
  117. enddo
  118. C
  119. CC - MODES DE CALCUL EN DEFORMATIONS "PLANES" GENERALISEES
  120. CC CAS 3D :
  121. IF (IDIM.EQ.3) THEN
  122. C RIEN FAIRE !
  123. CC CAS 2D :
  124. ELSE IF(IDIM.EQ.2) THEN
  125. CCCC CAS 2D AXISYMETRIQUE
  126. IF (IFOU.EQ.0) THEN
  127. R1=0.D0
  128. R2=0.D0
  129. DO 150 ID=1,NBNN
  130. R1=R1+SH1(1,ID)*XE1(1,ID)
  131. R2=R2+SH1(1,ID)*XE2(1,ID)
  132. 150 CONTINUE
  133. if (r1.lt.-xpetit/xzprec.or.r1.gt.xpetit/xzprec)then
  134. XJAC(3,3)=R2/R1
  135. else
  136. xjac(3,3)=xgrand*xzprec
  137. endif
  138. ELSE
  139. CCCC CAS 2D PLAN
  140. XJAC(3,3)=1.D0
  141. CCCC CAS 2D PLAN DEFO GENE
  142. * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(3,1) (cf. PIOCAP)
  143. IF (IFOU.EQ.-3) THEN
  144. XJAC(3,3) = XJAC(3,3) + XE2(3,1)
  145. ENDIF
  146. ENDIF
  147. CC CAS 1D :
  148. ELSE IF (IDIM.EQ.1) THEN
  149. CCCC CAS 1D PLAN
  150. IF (IFOU.GE.3 .AND. IFOU.LE.11) THEN
  151. XJAC(2,2) = 1.D0
  152. XJAC(3,3) = 1.D0
  153. * Rq : "Deplacement" UY du PTGENE est stocke dans XE2(2,1) (cf. PIOCAP)
  154. IF (IFOU.EQ.7 .OR. IFOU.EQ.8 .OR. IFOU.EQ.11) THEN
  155. XJAC(2,2) = XJAC(2,2) + XE2(2,1)
  156. ENDIF
  157. * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(3,1) (cf. PIOCAP)
  158. c* IF (IFOU.EQ.9 .OR. IFOU.EQ.10 .OR. IFOU.EQ.11) THEN
  159. IF (IFOU.GE.9) THEN
  160. XJAC(3,3) = XJAC(3,3) + XE2(3,1)
  161. ENDIF
  162. CCCC CAS 1D AXIS et SPHE
  163. ELSE
  164. CALL DISTRR(XE1,SH1,NBNN,R1)
  165. CALL DISTRR(XE2,SH1,NBNN,R2)
  166. if (r1.lt.-xpetit/xzprec.or.r1.gt.xpetit/xzprec)then
  167. FR2R1 = R2 /R1
  168. else
  169. FR2R1=XGRAND*XZPREC
  170. endif
  171. XJAC(3,3) = FR2R1
  172. CCCC CAS 1D SPHE
  173. IF (IFOU.EQ.15) THEN
  174. XJAC(2,2) = FR2R1
  175. ELSE
  176. CCCC CAS 1D AXIS
  177. XJAC(2,2) = 1.D0
  178. * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(2,1) (cf. PIOCAP)
  179. IF (IFOU.EQ.14) THEN
  180. XJAC(2,2) = XJAC(2,2) + XE2(2,1)
  181. ENDIF
  182. ENDIF
  183. ENDIF
  184. ENDIF
  185. C
  186. GO TO (500,600,700),KCAS
  187. C
  188. C KCAS=1 CAS DES CONTRAINTES
  189. C ----------------------------
  190. C
  191. 500 CONTINUE
  192. C
  193. CCCCCCCCCCCC CALCUL DE L'INVERSE DU DETERMINANT DE F
  194. C
  195. IF (IDIM.EQ.3) THEN
  196. DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  197. DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  198. DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  199. ELSE IF (IDIM.EQ.2) THEN
  200. DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3)
  201. ELSE IF (IDIM.EQ.1) THEN
  202. DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3)
  203. ENDIF
  204. if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  205. DETF=1./(DETF)
  206. else
  207. DETF=XGRAND*xzprec
  208. endif
  209. C
  210. C CALCUL DES CONTRAINTES DE CAUCHY
  211. C
  212. DO 160 ID=1,NCOELE
  213. IND=IN(ID)
  214. JND=JN(ID)
  215. TABCD = TAB(IC,ID)
  216. DO 170 IE=1,IDIM
  217. TABAU1=0.D0
  218. DO 171 IF=1,IDIM
  219. ICO=ITAB(IF,IE)
  220. TABAU1=TAB1(IC,ICO)*XJAC(JND,IF)
  221. 1 +TABAU1
  222. 171 CONTINUE
  223. TABAU1 = TABAU1 * XJAC(IND,IE)
  224. TABCD = TABCD + TABAU1
  225. 170 CONTINUE
  226. TAB(IC,ID)=TABCD*DETF + TAB(IC,ID)
  227. 160 CONTINUE
  228. C
  229. C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  230. C
  231. IF (IDIM.EQ.2) THEN
  232. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  233. ELSE IF (IDIM.EQ.1) THEN
  234. TAB(IC,2)=TAB1(IC,2)*XJAC(2,2)*XJAC(2,2)*DETF
  235. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  236. ENDIF
  237. GO TO 130
  238. C
  239. C KCAS=2 CAS DES DEFORMATIONS
  240. C -----------------------------
  241. C
  242. 600 CONTINUE
  243. C
  244. C
  245. CCCCCCCCCCCC CALCUL DE L'INVERSE DE F
  246. C
  247. CALL INVMA1(XJAC,3,3,KERRE)
  248. IF(KERRE.NE.0) THEN
  249. WRITE(6,77881) ((XJAC(MI,MJ),MJ=1,3),MI=1,3)
  250. 77881 FORMAT(2X,' MATRICE SINGULIERE' /(3(1X,1PE12.5)/))
  251. RETURN
  252. ENDIF
  253. C
  254. C CALCUL DES DEFORMATIONS
  255. C
  256. DO 260 ID=1,NCOELE
  257. IND=IN(ID)
  258. JND=JN(ID)
  259. TABCD=0.D0
  260. DO 271 IE=1,IDIM
  261. TABAU1=0.D0
  262. DO 270 IF=1,IDIM
  263. ICO=ITAB(IF,IE)
  264. TABAU1=TABAU1+
  265. 1 FAC(ICO)*TAB1(IC,ICO)*XJAC(IF,JND)
  266. 270 CONTINUE
  267. TABCD=TABCD+TABAU1*XJAC(IE,IND)
  268. 271 CONTINUE
  269. TAB(IC,ID)=TAB(IC,ID)+TABCD*FACI(ID)
  270. 260 CONTINUE
  271. C
  272. C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  273. C
  274. IF(IDIM.EQ.2) THEN
  275. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)
  276. ELSE IF(IDIM.EQ.1) THEN
  277. TAB(IC,2)=TAB1(IC,2)*XJAC(2,2)*XJAC(2,2)*DETF
  278. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  279. ENDIF
  280. C
  281. GO TO 130
  282. C
  283. C KCAS=3 CAS DE LA MATRICE DE HOOKE
  284. C -----------------------------------
  285. C
  286. 700 CONTINUE
  287. C
  288. CCCCCCCCCCCC CALCUL DE L'INVERSE DU DETERMINANT DE F
  289. C
  290. IF (IDIM.EQ.3) THEN
  291. DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  292. DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  293. DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  294. ELSE IF (IDIM.EQ.2) THEN
  295. DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3)
  296. ELSE IF (IDIM.EQ.1) THEN
  297. DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3)
  298. ENDIF
  299. if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  300. DETF=1./(DETF)
  301. else
  302. DETF=XGRAND*xzprec
  303. endif
  304. C
  305. IJ=1
  306. DO JJ=1,LHOOK
  307. DO II=1,LHOOK
  308. DDHOOK(II,JJ)=PRODDI(IC,IJ)
  309. IJ=IJ+1
  310. enddo
  311. enddo
  312. *
  313. CALL ZERO(DDHOMU,LHOOK,LHOOK)
  314. KEY=2
  315. C
  316. DO 720 LC=1,LHOOK
  317.  
  318. CALL ZERO (VEC,LHOOK,1)
  319. DO 760 ID=1,LHOOK
  320. IND=IN(ID)
  321. JND=JN(ID)
  322. DO IE=1,IDIM
  323. DO IF=1,IDIM
  324. ICO=ITAB(IF,IE)
  325. IF(ICO.EQ.LC) THEN
  326. VEC(ID)=VEC(ID)+
  327. & FAC(ICO)*XJAC(IE,IND)*XJAC(IF,JND)
  328. ENDIF
  329. enddo
  330. enddo
  331. VEC(ID)=VEC(ID)*FACI(ID)
  332. 760 CONTINUE
  333. C
  334. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  335. C
  336. IF (IDIM.EQ.2) THEN
  337. IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3)
  338. *
  339. ELSE IF (IDIM.EQ.1) THEN
  340. IF(LC.EQ.2) VEC(2)=XJAC(2,2)*XJAC(2,2)
  341. IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3)
  342. ENDIF
  343. C
  344. CALL MATVE1(DDHOOK,VEC,LHOOK,LHOOK,VEC2,KEY)
  345. C
  346. DO 761 ID=1,LHOOK
  347. IND=IN(ID)
  348. JND=JN(ID)
  349. DO IE=1,IDIM
  350. DO IF=1,IDIM
  351. ICO=ITAB(IF,IE)
  352. DDHOMU(ID,LC)=VEC2(ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF
  353. 1 +DDHOMU(ID,LC)
  354. enddo
  355. enddo
  356. 761 CONTINUE
  357. C
  358. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  359. C
  360. IF (IDIM.EQ.2) THEN
  361. DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  362. ELSE IF (IDIM.EQ.1) THEN
  363. DDHOMU(2,LC)=VEC2(2)*XJAC(2,2)*XJAC(2,2)*DETF
  364. DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  365. ENDIF
  366.  
  367. 720 CONTINUE
  368. C
  369. IJ=1
  370. DO JJ=1,LHOOK
  371. DO II=1,LHOOK
  372. PRODDO(IC,IJ)=DDHOMU(II,JJ)
  373. IJ=IJ+1
  374. enddo
  375. enddo
  376. *
  377. *
  378. 130 CONTINUE
  379. RETURN
  380. END
  381.  
  382.  
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  

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