Télécharger piocaf.eso

Retour à la liste

Numérotation des lignes :

piocaf
  1. C PIOCAF SOURCE PV090527 24/03/14 21:15:02 11870
  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. * unrolling de la boucle 270 pour mettre xjac en registre
  262. ** TABAU1=0.D0
  263. ** DO 270 IF=1,IDIM
  264. ** ICO=ITAB(IF,IE)
  265. ** TABAU1=TABAU1+
  266. ** 1 FAC(ICO)*TAB1(IC,ICO)*XJAC(IF,JND)
  267. *270 CONTINUE
  268. ICO1=ITAB(1,IE)
  269. tabau1= FAC(ICO1)*TAB1(IC,ICO1)*XJAC(1,JND)
  270. ICO2=ITAB(2,IE)
  271. tabau1= tabau1 + FAC(ICO2)*TAB1(IC,ICO2)*XJAC(2,JND)
  272. if(idim.eq.3) then
  273. ICO3=ITAB(3,IE)
  274. tabau1= tabau1 + FAC(ICO3)*TAB1(IC,ICO3)*XJAC(3,JND)
  275. endif
  276. TABCD=TABCD+TABAU1*XJAC(IE,IND)
  277. 271 CONTINUE
  278. TAB(IC,ID)=TAB(IC,ID)+TABCD*FACI(ID)
  279. 260 CONTINUE
  280. C
  281. C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  282. C
  283. IF(IDIM.EQ.2) THEN
  284. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)
  285. ELSE IF(IDIM.EQ.1) THEN
  286. TAB(IC,2)=TAB1(IC,2)*XJAC(2,2)*XJAC(2,2)*DETF
  287. TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  288. ENDIF
  289. C
  290. GO TO 130
  291. C
  292. C KCAS=3 CAS DE LA MATRICE DE HOOKE
  293. C -----------------------------------
  294. C
  295. 700 CONTINUE
  296. C
  297. CCCCCCCCCCCC CALCUL DE L'INVERSE DU DETERMINANT DE F
  298. C
  299. IF (IDIM.EQ.3) THEN
  300. DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  301. DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  302. DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  303. ELSE IF (IDIM.EQ.2) THEN
  304. DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3)
  305. ELSE IF (IDIM.EQ.1) THEN
  306. DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3)
  307. ENDIF
  308. if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  309. DETF=1./(DETF)
  310. else
  311. DETF=XGRAND*xzprec
  312. endif
  313. C
  314. IJ=1
  315. DO JJ=1,LHOOK
  316. DO II=1,LHOOK
  317. DDHOOK(II,JJ)=PRODDI(IC,IJ)
  318. IJ=IJ+1
  319. enddo
  320. enddo
  321. *
  322. CALL ZERO(DDHOMU,LHOOK,LHOOK)
  323. KEY=2
  324. C
  325. DO 720 LC=1,LHOOK
  326.  
  327. CALL ZERO (VEC,LHOOK,1)
  328. DO 760 ID=1,LHOOK
  329. IND=IN(ID)
  330. JND=JN(ID)
  331. DO IE=1,IDIM
  332. DO IF=1,IDIM
  333. ICO=ITAB(IF,IE)
  334. IF(ICO.EQ.LC) THEN
  335. VEC(ID)=VEC(ID)+
  336. & FAC(ICO)*XJAC(IE,IND)*XJAC(IF,JND)
  337. ENDIF
  338. enddo
  339. enddo
  340. VEC(ID)=VEC(ID)*FACI(ID)
  341. 760 CONTINUE
  342. C
  343. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  344. C
  345. IF (IDIM.EQ.2) THEN
  346. IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3)
  347. *
  348. ELSE IF (IDIM.EQ.1) THEN
  349. IF(LC.EQ.2) VEC(2)=XJAC(2,2)*XJAC(2,2)
  350. IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3)
  351. ENDIF
  352. C
  353. CALL MATVE1(DDHOOK,VEC,LHOOK,LHOOK,VEC2,KEY)
  354. C
  355. DO 761 ID=1,LHOOK
  356. IND=IN(ID)
  357. JND=JN(ID)
  358. DO IE=1,IDIM
  359. DO IF=1,IDIM
  360. ICO=ITAB(IF,IE)
  361. DDHOMU(ID,LC)=VEC2(ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF
  362. 1 +DDHOMU(ID,LC)
  363. enddo
  364. enddo
  365. 761 CONTINUE
  366. C
  367. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  368. C
  369. IF (IDIM.EQ.2) THEN
  370. DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  371. ELSE IF (IDIM.EQ.1) THEN
  372. DDHOMU(2,LC)=VEC2(2)*XJAC(2,2)*XJAC(2,2)*DETF
  373. DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  374. ENDIF
  375.  
  376. 720 CONTINUE
  377. C
  378. IJ=1
  379. DO JJ=1,LHOOK
  380. DO II=1,LHOOK
  381. PRODDO(IC,IJ)=DDHOMU(II,JJ)
  382. IJ=IJ+1
  383. enddo
  384. enddo
  385. *
  386. *
  387. 130 CONTINUE
  388. RETURN
  389. END
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  

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