Télécharger piocaf.eso

Retour à la liste

Numérotation des lignes :

  1. C PIOCAF SOURCE PV 18/11/21 21:15:08 10002
  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)
  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 IN/1,2,3,1,1,2/
  74. DATA JN/1,2,3,2,3,3/
  75. C
  76. DATA ITAB(1,1),ITAB(1,2),ITAB(1,3)/1,4,5/
  77. DATA ITAB(2,1),ITAB(2,2),ITAB(2,3)/4,2,6/
  78. DATA ITAB(3,1),ITAB(3,2),ITAB(3,3)/5,6,3/
  79. C
  80. KERRE=0
  81. C
  82. C MISE A ZERO DES CONTRAINTES OU DES DEFORMATIONS
  83. C
  84. ** DO 50 IA=1,NBPTEL
  85. ** DO 50 IB=1,NCOELE
  86. ** TAB(IA,IB)=0.D0
  87. **50 CONTINUE
  88. CALL ZERO(TAB,NBPTEL,NCOELE)
  89. C
  90. C BOUCLE SUR LES POINTS DE GAUSS
  91. C
  92. DO 130 IC=1,NBPTEL
  93. C
  94. if (nbsh.ne.nbnn) then
  95. CALL HPRIMEX(XE1,NBNN,IDIM,NBSH,SHP,IC,SH1,DJAC)
  96. else
  97. CALL HPRIME(XE1,NBNN,IDIM,SHP,IC,SH1,DJAC)
  98. endif
  99. C
  100. C CALCUL DE LA MATRICE F (=XJAC)
  101. C
  102. CALL ZERO(XJAC,3,3)
  103. DO IF=1,IDIM
  104. JF = IF + 1
  105. r_z1 = 0.D0
  106. r_z2 = 0.D0
  107. r_z3 = 0.D0
  108. do id=1,nbnn
  109. r_z1 = r_z1 + SH1(JF,ID)*XE2(1,ID)
  110. if (idim.ge.2) r_z2 = r_z2 + SH1(JF,ID)*XE2(2,ID)
  111. if (idim.ge.3) r_z3 = r_z3 + SH1(JF,ID)*XE2(3,ID)
  112. enddo
  113. if (idim.ge.1) XJAC(1,IF) = r_z1
  114. if (idim.ge.2) XJAC(2,IF) = r_z2
  115. if (idim.ge.3) XJAC(3,IF) = r_z3
  116. enddo
  117. C
  118. CC - MODES DE CALCUL EN DEFORMATIONS "PLANES" GENERALISEES
  119. CC CAS 3D :
  120. IF (IDIM.EQ.3) THEN
  121. C RIEN FAIRE !
  122. CC CAS 2D :
  123. ELSE IF(IDIM.EQ.2) THEN
  124. CCCC CAS 2D AXISYMETRIQUE
  125. IF (IFOU.EQ.0) THEN
  126. R1=0.D0
  127. R2=0.D0
  128. DO 150 ID=1,NBNN
  129. R1=R1+SH1(1,ID)*XE1(1,ID)
  130. R2=R2+SH1(1,ID)*XE2(1,ID)
  131. 150 CONTINUE
  132. if (r1.lt.-xpetit/xzprec.or.r1.gt.xpetit/xzprec)then
  133. XJAC(3,3)=R2/R1
  134. else
  135. xjac(3,3)=xgrand*xzprec
  136. endif
  137. ELSE
  138. CCCC CAS 2D PLAN
  139. XJAC(3,3)=1.D0
  140. CCCC CAS 2D PLAN DEFO GENE
  141. * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(3,1) (cf. PIOCAP)
  142. IF (IFOU.EQ.-3) THEN
  143. XJAC(3,3) = XJAC(3,3) + XE2(3,1)
  144. ENDIF
  145. ENDIF
  146. CC CAS 1D :
  147. ELSE IF (IDIM.EQ.1) THEN
  148. CCCC CAS 1D PLAN
  149. IF (IFOU.GE.3 .AND. IFOU.LE.11) THEN
  150. XJAC(2,2) = 1.D0
  151. XJAC(3,3) = 1.D0
  152. * Rq : "Deplacement" UY du PTGENE est stocke dans XE2(2,1) (cf. PIOCAP)
  153. IF (IFOU.EQ.7 .OR. IFOU.EQ.8 .OR. IFOU.EQ.11) THEN
  154. XJAC(2,2) = XJAC(2,2) + XE2(2,1)
  155. ENDIF
  156. * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(3,1) (cf. PIOCAP)
  157. c* IF (IFOU.EQ.9 .OR. IFOU.EQ.10 .OR. IFOU.EQ.11) THEN
  158. IF (IFOU.GE.9) THEN
  159. XJAC(3,3) = XJAC(3,3) + XE2(3,1)
  160. ENDIF
  161. CCCC CAS 1D AXIS et SPHE
  162. ELSE
  163. CALL DISTRR(XE1,SH1,NBNN,R1)
  164. CALL DISTRR(XE2,SH1,NBNN,R2)
  165. if (r1.lt.-xpetit/xzprec.or.r1.gt.xpetit/xzprec)then
  166. FR2R1 = R2 /R1
  167. else
  168. FR2R1=XGRAND*XZPREC
  169. endif
  170. XJAC(3,3) = FR2R1
  171. CCCC CAS 1D SPHE
  172. IF (IFOU.EQ.15) THEN
  173. XJAC(2,2) = FR2R1
  174. ELSE
  175. CCCC CAS 1D AXIS
  176. XJAC(2,2) = 1.D0
  177. * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(2,1) (cf. PIOCAP)
  178. IF (IFOU.EQ.14) THEN
  179. XJAC(2,2) = XJAC(2,2) + XE2(2,1)
  180. ENDIF
  181. ENDIF
  182. ENDIF
  183. ENDIF
  184. C
  185. GO TO (500,600,700),KCAS
  186. C
  187. C KCAS=1 CAS DES CONTRAINTES
  188. C ----------------------------
  189. C
  190. 500 CONTINUE
  191. C
  192. CCCCCCCCCCCC CALCUL DE L'INVERSE DU DETERMINANT DE F
  193. C
  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. ELSE IF (IDIM.EQ.2) THEN
  199. DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3)
  200. ELSE IF (IDIM.EQ.1) THEN
  201. DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3)
  202. ENDIF
  203. if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  204. DETF=1./(DETF)
  205. else
  206. DETF=XGRAND*xzprec
  207. endif
  208. C
  209. C CALCUL DES CONTRAINTES DE CAUCHY
  210. C
  211. DO 160 ID=1,NCOELE
  212. IND=IN(ID)
  213. JND=JN(ID)
  214. TABCD = TAB(IC,ID)
  215. DO 170 IE=1,IDIM
  216. TABAU1=0.D0
  217. DO 171 IF=1,IDIM
  218. ICO=ITAB(IE,IF)
  219. TABAU1=TAB1(IC,ICO)*XJAC(JND,IF)
  220. 1 +TABAU1
  221. 171 CONTINUE
  222. TABAU1 = TABAU1 * XJAC(IND,IE)
  223. TABCD = TABCD + TABAU1
  224. 170 CONTINUE
  225. TAB(IC,ID)=TABCD*DETF + TAB(IC,ID)
  226. 160 CONTINUE
  227. C
  228. C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  229. C
  230. IF (IDIM.EQ.2) THEN
  231. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  232. ELSE IF (IDIM.EQ.1) THEN
  233. TAB(IC,2)=TAB1(IC,2)*XJAC(2,2)*XJAC(2,2)*DETF
  234. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  235. ENDIF
  236. GO TO 130
  237. C
  238. C KCAS=2 CAS DES DEFORMATIONS
  239. C -----------------------------
  240. C
  241. 600 CONTINUE
  242. C
  243. C
  244. CCCCCCCCCCCC CALCUL DE L'INVERSE DE F
  245. C
  246. CALL INVMA1(XJAC,3,3,KERRE)
  247. IF(KERRE.NE.0) THEN
  248. WRITE(6,77881) ((XJAC(MI,MJ),MJ=1,3),MI=1,3)
  249. 77881 FORMAT(2X,' MATRICE SINGULIERE' /(3(1X,1PE12.5)/))
  250. RETURN
  251. ENDIF
  252. C
  253. C CALCUL DES DEFORMATIONS
  254. C
  255. DO 260 ID=1,NCOELE
  256. IND=IN(ID)
  257. JND=JN(ID)
  258. TABCD=0.D0
  259. DO 271 IE=1,IDIM
  260. TABAU1=0.D0
  261. DO 270 IF=1,IDIM
  262. ICO=ITAB(IE,IF)
  263. TABAU1=TABAU1+
  264. 1 FAC(ICO)*TAB1(IC,ICO)*XJAC(IF,JND)
  265. 270 CONTINUE
  266. TABCD=TABCD+TABAU1*XJAC(IE,IND)
  267. 271 CONTINUE
  268. TAB(IC,ID)=TAB(IC,ID)+TABCD/FAC(ID)
  269. 260 CONTINUE
  270. C
  271. C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  272. C
  273. IF(IDIM.EQ.2) THEN
  274. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)
  275. ELSE IF(IDIM.EQ.1) THEN
  276. TAB(IC,2)=TAB1(IC,2)*XJAC(2,2)*XJAC(2,2)*DETF
  277. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  278. ENDIF
  279. C
  280. GO TO 130
  281. C
  282. C KCAS=3 CAS DE LA MATRICE DE HOOKE
  283. C -----------------------------------
  284. C
  285. 700 CONTINUE
  286. C
  287. CCCCCCCCCCCC CALCUL DE L'INVERSE DU DETERMINANT DE F
  288. C
  289. IF (IDIM.EQ.3) THEN
  290. DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  291. DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  292. DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  293. ELSE IF (IDIM.EQ.2) THEN
  294. DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3)
  295. ELSE IF (IDIM.EQ.1) THEN
  296. DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3)
  297. ENDIF
  298. if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  299. DETF=1./(DETF)
  300. else
  301. DETF=XGRAND*xzprec
  302. endif
  303. C
  304. IJ=1
  305. DO 710 JJ=1,LHOOK
  306. DO 710 II=1,LHOOK
  307. DDHOOK(II,JJ)=PRODDI(IC,IJ)
  308. IJ=IJ+1
  309. 710 CONTINUE
  310. *
  311. CALL ZERO(DDHOMU,LHOOK,LHOOK)
  312. KEY=2
  313. C
  314. DO 720 LC=1,LHOOK
  315.  
  316. CALL ZERO (VEC,LHOOK,1)
  317. DO 760 ID=1,LHOOK
  318. IND=IN(ID)
  319. JND=JN(ID)
  320. DO 770 IE=1,IDIM
  321. DO 770 IF=1,IDIM
  322. ICO=ITAB(IE,IF)
  323. IF(ICO.EQ.LC) THEN
  324. VEC(ID)=VEC(ID)+
  325. & FAC(ICO)*XJAC(IE,IND)*XJAC(IF,JND)/FAC(ID)
  326. ENDIF
  327. 770 CONTINUE
  328. 760 CONTINUE
  329. C
  330. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  331. C
  332. IF (IDIM.EQ.2) THEN
  333. IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3)
  334. *
  335. ELSE IF (IDIM.EQ.1) THEN
  336. IF(LC.EQ.2) VEC(2)=XJAC(2,2)*XJAC(2,2)
  337. IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3)
  338. ENDIF
  339. C
  340. CALL MATVE1(DDHOOK,VEC,LHOOK,LHOOK,VEC2,KEY)
  341. C
  342. DO 761 ID=1,LHOOK
  343. IND=IN(ID)
  344. JND=JN(ID)
  345. DO 771 IE=1,IDIM
  346. DO 771 IF=1,IDIM
  347. ICO=ITAB(IE,IF)
  348. DDHOMU(ID,LC)=VEC2(ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF
  349. 1 +DDHOMU(ID,LC)
  350. 771 CONTINUE
  351. 761 CONTINUE
  352. C
  353. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  354. C
  355. IF (IDIM.EQ.2) THEN
  356. DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  357. ELSE IF (IDIM.EQ.1) THEN
  358. DDHOMU(2,LC)=VEC2(2)*XJAC(2,2)*XJAC(2,2)*DETF
  359. DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  360. ENDIF
  361.  
  362. 720 CONTINUE
  363. C
  364. IJ=1
  365. DO 780 JJ=1,LHOOK
  366. DO 780 II=1,LHOOK
  367. PRODDO(IC,IJ)=DDHOMU(II,JJ)
  368. IJ=IJ+1
  369. 780 CONTINUE
  370. *
  371. *
  372. 130 CONTINUE
  373. RETURN
  374. END
  375.  
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  

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