Télécharger la plaquette Cast3M


fr en

Télécharger shb8.eso

Retour à la liste

Numérotation des lignes : 

shb8

  1. C SHB8      SOURCE    PV        18/06/18    21:15:29     9860           
  2.       SUBROUTINE SHB8(ICLE,XE,D,PROPEL,WORK,RE,OUT)
  3.       IMPLICIT REAL *8(A-H,O-Z)
  4.  
  5. *      SUBROUTINE SHB8(ICLE,XE,D,PROPEL,WORK,RE,OUT)
  6.  
  7. C
  8. C               ELEMENT SHB8
  9. C
  10.  
  11.       INTEGER IOP,IBID,IDPLA(5),LAG,IRDC
  12.       DIMENSION EYG(5),ENU(5)
  13.       DIMENSION CMATTMP(6,6)
  14.       DIMENSION XE(24),D(*),PROPEL(*),WORK(30),RE(24,24),OUT(*)
  15.       DIMENSION XXG5(5),PXG5(5),XCOQ(3,4),BKSIP(3,8,5),B(3,8)
  16.       DIMENSION XELOC(24),XCENT(3),PPP(3,3),PPPT(3,3)
  17.       DIMENSION XL(3,4),XXX(3),YYY(3),BGL_LOC(6,24)
  18.       DIMENSION BGL_LOCT(24,6),TMPTAB(6,24),TMPKE(24,24),CMATLOC(6,6)
  19.       DIMENSION XELOCTMP(24)
  20.       DIMENSION XXVB(3),HIJ(6),XKSTAB(24,24),TMPKE2(24,24)
  21.       DIMENSION GB(8,4),GS(8,4),VB(8,3),XXGB(3,4)
  22.       DIMENSION XK11(8,8),XK22(8,8),XK33(8,8),RR2(3,3),XK12(8,8)
  23.       DIMENSION XK21(8,8),XK13(8,8),XK23(8,8),XK31(8,8),XK32(8,8)
  24.       DIMENSION XR(3,3),XRT(3,3)
  25. C
  26.       DIMENSION DEPS(6),DUSDX(9),UE(3,8),SIG(6),XE_LOCP(24)
  27.       DIMENSION UDEF(24),XE_LOC(24),XE_LOC12(24)
  28.       DIMENSION DEPS_LOC(6),SIG_LOC(6)
  29. C
  30.       DIMENSION F(3,8),BSIG(24),SIGMAG(6),ULOC(3,8),PQIALPH(3,4)
  31.       DIMENSION QIALPHA(3,4),FHG(3,8),RR12(3,3),RR1(3,3),FHG24(24)
  32. C
  33.       DIMENSION DIREC(3),V(3),X56(3),X67(3),X78(3),X85(3),X12(3)
  34.       DIMENSION X23(3),X34(3),X41(3),V56_78(3),V85_67(3),VSUP(3)
  35.       DIMENSION V12_34(3),V41_23(3),VINF(3)
  36. C
  37.       DIMENSION XKSIGTMP1(8,8),XKSIGTMP2(8,8)
  38. C
  39.       DIMENSION XCOQ14(3),XCOQ23(3),XCOQ34(3),XCOQ12(3),XNORMALE(3)
  40.       DIMENSION XCOQ12_34(3),XCOQ14_23(3),XKP(24,24),BKP(3,4)
  41. C
  42.       DIMENSION XMASSE(8,8)
  43. C
  44.       DIMENSION FQ(24)
  45. C
  46. C     AVRIL 2004: PROGRAMMATION DU SHB8 AVEC B_2_CHAPEAU_BAR
  47. C
  48.       DIMENSION CB2BAR(6,6),SIGMOY(6)
  49. C
  50. CCCCCCCCCCCCCC ENTREES CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  51. C          ICLE=2    ON CALCULE LA MATRICE DE RAIDEUR
  52. C          ICLE=3    ON CALCULE LA MATRICE DE MASSE
  53. C          ICLE=4    ON CALCULE LA MATRICE TANGENTE
  54. C          ICLE=5    ON CALCULE LES FORCES DE PRESSION
  55. C          ICLE=7    ON CALCULE LES CONTRAINTES
  56. C          ICLE=8    ON CALCULE BSIGMA
  57. C          ICLE=9    ON CALCULE KSIGMA
  58. C          ICLE=10   ON CALCULE KP
  59. C          ICLE=11   On calcule epsi a partir de deplacement
  60. C          icle=12   on calcule epsi à partir des contraintes
  61. C                    ou  sigma à partir de epsi suivan propel(3)
  62. C
  63. C          D         MATRICE D'ELASTICITE
  64. C          PROPEL    PROPRIETES ELASTIQUES (1)=YUNG (2)=NU
  65. C                                          (3)=RO   (4)=ALPHA
  66. C          WORK      TABLEAU DE TRAVAIL
  67. CCCCCCCCCCCCCC SORTIE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  68. C
  69. C          RE        MATRICE
  70. C          OUT       VECTEUR FORCE OU VECTEUR CONTRAINTE AUX NOEUDS
  71. C
  72. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  73. C
  74. C INITIALISATIONS
  75. C
  76. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  77. C INFOS:
  78. C XE EST RANGE COMME CA: (XNOEUD1 YNOEUD1 ZNOEUD1, XNOEUD2 YNOEUD2 ZNOEUD2,...)
  79. C DANS SHB8_TEST_NUM: ATTENTION A LA NUMEROTATION DES NOEUDS
  80. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  81.       DATA GB/1.D0, 1.D0,-1.D0,-1.D0,-1.D0,-1.D0, 1.D0, 1.D0,
  82.      &        1.D0,-1.D0,-1.D0, 1.D0,-1.D0, 1.D0, 1.D0,-1.D0,
  83.      &        1.D0,-1.D0, 1.D0,-1.D0, 1.D0,-1.D0, 1.D0,-1.D0,
  84.      &       -1.D0, 1.D0,-1.D0, 1.D0, 1.D0,-1.D0, 1.D0,-1.D0/
  85. C
  86. C VB: COORD DES NOEUDS DANS REPERE DE REFERENCE
  87. C
  88.       DATA VB/-1.D0, 1.D0, 1.D0,-1.D0,-1.D0, 1.D0, 1.D0,-1.D0,
  89.      &        -1.D0,-1.D0, 1.D0, 1.D0,-1.D0,-1.D0, 1.D0, 1.D0,
  90.      &        -1.D0,-1.D0,-1.D0,-1.D0, 1.D0, 1.D0, 1.D0, 1.D0/
  91. C
  92. C ON DEFINI LES POINTS GAUSS ET LES POIDS
  93. C
  94.       XXG5(1) = -0.906179845938664D0
  95.       XXG5(2) = -0.538469310105683D0
  96.       XXG5(3) =  0.D0
  97.       XXG5(4) =  0.538469310105683D0
  98.       XXG5(5) =  0.906179845938664D0
  99. C
  100.       PXG5(1) =  0.236926885056189D0
  101.       PXG5(2) =  0.478628670499366D0
  102.       PXG5(3) =  0.568888888888889D0
  103.       PXG5(4) =  0.478628670499366D0
  104.       PXG5(5) =  0.236926885056189D0
  105. C
  106.       UNS8    = 1.D0/8.D0
  107.       UNS3    = 1.D0/3.D0
  108.       UNS5    = 1.D0/5.D0
  109.       UNS6    = 1.D0/6.D0
  110.       ZERO    = 0.D0
  111. C
  112. C TYPE DE LOI DE COMPORTEMENT:
  113. C     IRDC = 1 : SHB8 TYPE PLEXUS
  114. C     IRDC = 2 : C.P.
  115. C     IRDC = 3 : 3D COMPLETE
  116. C
  117. C      IRDC = 1
  118.       IRDC=nint(OUT(1))
  119.       IF(ICLE.EQ.2)THEN
  120. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  121. C                                                                  C
  122. C ON CALCULE LA RAIDEUR : SORTIE DANS RE                           C
  123. C                                                                  C
  124. C SI IETAN = 1 , ALORS ON CALCULE AUSSI                            C
  125. C                LA MATRICE TANGENTE PLASTIQUE                     C
  126. C                                                                  C
  127. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  128.       IETAN   = nint(PROPEL(13))
  129. C
  130. C INTIALISATION LONGUEUR DES COTES
  131. C CALCUL DES COEFF D ELANCEMENT A METTRE DANS LA MATRICE DE CPT
  132. C
  133.       XXL1 = 0.D0
  134.       XXL2 = 0.D0
  135.       TT1  = 0.D0
  136.       TT2  = 0.D0
  137. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  138. C
  139. C STABILISATION ADAPTATIVE EN FONCTION DE LA DISTORTION DE L'ELEMENT
  140. C
  141. Cbp  POUR L'INSTANT, ON NE MET PAS EN SERVICE => on commente
  142. Cbp       DO I=1,3
  143. Cbp C DISTANCE ENTRE 1 ET 5 (EPAISSEUR)
  144. Cbp          TT1   = TT1+(XE(I+12)-XE(I))**2
  145. Cbp C DISTANCE ENTRE 3 ET 7 (EPAISSEUR)
  146. Cbp          TT2   = TT2+(XE(I+18)-XE(I+6))**2
  147. Cbp C DISTANCE ENTRE 1 ET 2
  148. Cbp          XXL1  = XXL1+(XE(I+3)-XE(I))**2
  149. Cbp C DISTANCE ENTRE 2 ET 3
  150. Cbp          XXL2  = XXL2+(XE(I+6)-XE(I+3))**2
  151. Cbp       ENDDO
  152. Cbp       XXL1   = SQRT(XXL1)
  153. Cbp       XXL2   = SQRT(XXL2)
  154. Cbp       TT1    = 0.5*(SQRT(TT1)+SQRT(TT2))
  155. Cbp       COELA1 = 5./6.
  156. Cbp       COELA2 = 5./6.
  157. Cbp C ELANCEMENT DANS DIRECTION 2
  158. Cbp       ELT    = 6.*TT1/XXL1
  159. Cbp       IF(COELA1.GT.ELT) COELA1=ELT
  160. Cbp C ELANCEMENT DANS DIRECTION 1
  161. Cbp       ELT    = 6.*TT1/XXL2
  162. Cbp       IF(COELA2.GT.ELT) COELA2=ELT
  163. C POUR L'INSTANT, ON NE MET PAS EN SERVICE:
  164.       COELA1 = 1.D0
  165.       COELA2 = 1.D0
  166. C
  167. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  168. C
  169.       IF (IETAN.EQ.1) THEN
  170.          EYG(1)  = PROPEL(3)
  171.          EYG(2)  = PROPEL(4)
  172.          EYG(3)  = PROPEL(5)
  173.          EYG(4)  = PROPEL(6)
  174.          EYG(5)  = PROPEL(7)
  175.          ENU(1)  = PROPEL(8)
  176.          ENU(2)  = PROPEL(9)
  177.          ENU(3)  = PROPEL(10)
  178.          ENU(4)  = PROPEL(11)
  179.          ENU(5)  = PROPEL(12)
  180.          DO I=1,5
  181. C     SI IDPLA(I)=1, POINT GAUSS (I) PLASTIQUE
  182.             IDPLA(I)=nint(RE(I,1))
  183.          ENDDO
  184.       ENDIF
  185. C
  186.       CALL ZDANUL(RE,576)
  187.       CALL ZDANUL(CMATLOC,36)
  188. C      CALL ZDANUL(CMATTMP,36)
  189.       IELEM=nint(WORK(1))
  190. C
  191. C PROCEDURE DE TEST DE NUMEROTATION DES ELEMENTS: MARCHE BIEN, MAIS DEJA FAIT EN
  192. C                                                 AMONT
  193. C
  194. C       CALL SHB8_TEST_NUM(XE,IELEM)
  195. C
  196. C ON DEFINI CMATLOC: MATRICE DE COMPORTEMENT
  197. C
  198.       IF (IRDC.LE.2) THEN
  199. C     COMPORTEMENT SHB8 PLEXUS ou C.P.
  200.         XNU          = PROPEL(2)
  201.         XNU_B        = XNU / (1-XNU)
  202.         XLAMBDA      = PROPEL(1)*PROPEL(2)/(1-PROPEL(2)*PROPEL(2))
  203.         XMU          = 0.5D0 * PROPEL(1) / (1+PROPEL(2))      
  204.         CMATLOC(1,1) = XLAMBDA+2*XMU
  205.         CMATLOC(1,2) = XLAMBDA
  206.         CMATLOC(2,1) = XLAMBDA
  207.         CMATLOC(2,2) = XLAMBDA+2*XMU
  208.         IF (IRDC.EQ.1) THEN
  209. C       COMPORTEMENT SHB8 PLEXUS
  210.            CMATLOC(3,3) = PROPEL(1)
  211.         ENDIF
  212.         IF (IRDC.EQ.2) THEN
  213. C       COMPORTEMENT C.P.
  214.            CMATLOC(3,3) = 0.D0
  215.         ENDIF
  216.         CMATLOC(4,4) = XMU
  217.         CMATLOC(5,5) = XMU
  218.         CMATLOC(6,6) = XMU
  219.       ENDIF
  220. C
  221.       IF (IRDC.EQ.3) THEN
  222. C COMPORTEMENT LOI TRIDIM MMC 3D
  223.         XNU          = PROPEL(2)
  224.         XCOOEF       = PROPEL(1)/((1+XNU)*(1-2*XNU))
  225.         CMATLOC(1,1) = (1-XNU)*XCOOEF
  226.         CMATLOC(2,2) = (1-XNU)*XCOOEF
  227.         CMATLOC(3,3) = (1-XNU)*XCOOEF
  228.         CMATLOC(1,2) = XNU*XCOOEF
  229.         CMATLOC(2,1) = XNU*XCOOEF
  230.         CMATLOC(1,3) = XNU*XCOOEF
  231.         CMATLOC(3,1) = XNU*XCOOEF
  232.         CMATLOC(2,3) = XNU*XCOOEF
  233.         CMATLOC(3,2) = XNU*XCOOEF
  234.         CMATLOC(4,4) = (1-2*XNU)*0.5*XCOOEF
  235.         CMATLOC(5,5) = (1-2*XNU)*0.5*XCOOEF
  236.         CMATLOC(6,6) = (1-2*XNU)*0.5*XCOOEF
  237.       ENDIF
  238.  
  239. Cbp   cas non prevu      
  240.       IF (IRDC.GT.3) CALL ERREUR(5)
  241. C
  242. C      CALL SHIFTD(CMATLOC,CMATTMP,36)
  243. C
  244. C CALCUL DE BKSIP(3,8,IP) DANS REPERE DE REFERENCE
  245. C      BKSIP(1,*,IP) = VECTEUR BX AU POINT GAUSS IP
  246. C      BKSIP(2,*,IP) = VECTEUR BY AU POINT GAUSS IP
  247. C      BKSIP(3,*,IP) = VECTEUR BZ AU POINT GAUSS IP
  248. C
  249.       CALL CUBKSI(5,XXG5,BKSIP)
  250. C
  251. C DEBUT DE LA BOUCLE SUR LES 5 PTS GAUSS
  252. C
  253.       DO IP=1,5
  254. C
  255. C SI POINT PLASTIQUE, ON ITERE AVEC LE MODULE TANGENT
  256. C
  257. C      IF (IETAN.EQ.1) THEN
  258. C      IF (IDPLA(IP).EQ.1) THEN
  259. C      CMATLOC(1,1) = CMATTMP(1,1)*EYG(IP)/PROPEL(1)
  260. C      CMATLOC(2,2) = CMATTMP(2,2)*EYG(IP)/PROPEL(1)
  261. C      CMATLOC(3,3) = CMATTMP(3,3)*EYG(IP)/PROPEL(1)
  262. C      CMATLOC(1,2) = CMATTMP(1,2)*EYG(IP)/PROPEL(1)
  263. C      CMATLOC(2,1) = CMATTMP(2,1)*EYG(IP)/PROPEL(1)
  264. C      CMATLOC(4,4) = CMATTMP(4,4)*EYG(IP)/PROPEL(1)
  265. C      CMATLOC(5,5) = CMATTMP(5,5)*EYG(IP)/PROPEL(1)
  266. C      CMATLOC(6,6) = CMATTMP(6,6)*EYG(IP)/PROPEL(1)
  267. C      ELSE
  268. C      CMATLOC(1,1) = CMATTMP(1,1)
  269. C      CMATLOC(2,2) = CMATTMP(2,2)
  270. C      CMATLOC(3,3) = CMATTMP(3,3)
  271. C      CMATLOC(1,2) = CMATTMP(1,2)
  272. C      CMATLOC(2,1) = CMATTMP(2,1)
  273. C      CMATLOC(4,4) = CMATTMP(4,4)
  274. C      CMATLOC(5,5) = CMATTMP(5,5)
  275. C      CMATLOC(6,6) = CMATTMP(6,6)
  276. C      ENDIF
  277. C      ENDIF
  278. C
  279. C DEFINITION DES 4 POINTS  COQUES
  280. C
  281.          ZETA  = XXG5(IP)
  282.          ZLAMB = 0.5*(1.-ZETA)
  283.          DO I=1,4
  284.             DO J=1,3
  285.                XCOQ(J,I) = ZLAMB*XE((I-1)*3+J)
  286.      &                   + (1.-ZLAMB)*XE((I-1+4)*3+J)
  287.             END DO
  288.          END DO
  289. C
  290. C CALCUL DE PPP 3*3 PASSAGE DE GLOBAL A LOCAL,COQUE
  291. C XCENT : COORD GLOBAL DU CENTRE DE L'ELEMENT
  292. C
  293.          CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,XBID)
  294. C
  295. C CALCUL DE B EN GLOBAL
  296. C
  297. C ATTENTION A L'ORDRE DE EPSILON:
  298. C  FARID DANS SON PAPIER: 11 22 33 12 13 23
  299. C  HARID DANS PLEXUS:     11 22 33 12 23 13
  300. C
  301. C
  302. C ON RAJOUTE LES TERMES H1,X . G1  , H2,X . G2
  303. C                   ET  H1,Y . G1  , H2,Y . G2
  304. C AVEC H1   = Y.Z    H2   = X.Z
  305. C DONC H1,X =0       H2,X = Z
  306. C ET   H1,Y =Z       H2,Y = 0
  307. C
  308. C DONC IL NE RESTE PLUS QU'A CALCULER G1 ET G2, ET A AJOUTER A BKSIP
  309. C
  310.          CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC)
  311. C
  312. C IL FAUT:  EPS_LOCAL=BGLOB .U_NODAL_GLOBAL
  313. C ON CALCULE BGL_LOC LA MATRICE B(6,24) UGLOBAL ---> EPS LOCAL
  314. C
  315.          CALL ASBGLB(BGL_LOC,B,PPP)
  316. C
  317. C IL FAUT TRANSPOSER BGL_LOC
  318. C
  319.          CALL AEQBT(BGL_LOCT,BGL_LOC,6,24)
  320. C
  321. C CALCUL DE LA MATRICE TANGENTE PLASTIQUE: VOIR FICHIER SHB8_MAT_TGE_PLAS.FFF
  322. C
  323. C
  324. C IL NE RESTE PLUS QU'A FAIRE: BGL_LOCT * C * BGL_LOC
  325. C
  326.          CALL ZDANUL(TMPTAB,144)
  327.          CALL ZDANUL(TMPKE,576)
  328.          CALL ZDANUL(TMPKE2,576)
  329. *         CALL MULMAT(6,6,24,CMATLOC,BGL_LOC,TMPTAB)
  330.          CALL MULMAT( TMPTAB,CMATLOC,BGL_LOC,6,24,6)
  331. *         CALL MULMAT(24,6,24,BGL_LOCT,TMPTAB,TMPKE2)
  332.          CALL MULMAT(TMPKE2,BGL_LOCT,TMPTAB,24,24,6)
  333. C
  334. C ASSEMBLAGE: KE=KE + POIDS*JACOBIAN*TMPKE
  335. C
  336.          DO J=1,8
  337.             DO I=1,24
  338.                TMPKE(I,(J-1)*3+1)=TMPKE2(I,J)
  339.                TMPKE(I,(J-1)*3+2)=TMPKE2(I,J+8)
  340.                TMPKE(I,(J-1)*3+3)=TMPKE2(I,J+16)
  341.             END DO
  342.          END DO
  343.          CALL ZDANUL(TMPKE2,576)
  344.          DO I=1,8
  345.             DO J=1,24
  346.                TMPKE2((I-1)*3+1,J)=TMPKE(I,J)
  347.                TMPKE2((I-1)*3+2,J)=TMPKE(I+8,J)
  348.                TMPKE2((I-1)*3+3,J)=TMPKE(I+16,J)
  349.             END DO
  350.          END DO
  351.          DO J=1,24
  352.             DO I=1,24
  353.                RE(I,J)=RE(I,J) + 4.*AJAC*PXG5(IP)*TMPKE2(I,J)
  354.             END DO
  355.          END DO
  356.       END DO
  357. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  358. C                                                                  C
  359. C MATRICE DE STABILISATION : PAS DE BOUCLE SUR LES POINTS DE GAUSS C
  360. C                                                                  C
  361. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  362. C
  363. C  ON DEFINIE LES CONSTANTS DE B_2_CHAPEAU_BAR
  364. C
  365.       CALL ZDANUL(CB2BAR,36)
  366.  
  367.       CB2BAR(1,1) = 1.D0
  368.       CB2BAR(1,2) = 1.D0
  369.       CB2BAR(1,3) = ZERO
  370.       CB2BAR(1,4) = ZERO
  371.       CB2BAR(1,6) = ZERO
  372. C
  373.       CB2BAR(2,1) = ZERO
  374.       CB2BAR(2,2) = ZERO
  375.       CB2BAR(2,3) = 1.D0
  376.       CB2BAR(2,4) = 1.D0
  377.       CB2BAR(2,6) = ZERO
  378. C
  379.       CB2BAR(3,1) = ZERO
  380.       CB2BAR(3,2) = ZERO
  381.       CB2BAR(3,3) = ZERO
  382.       CB2BAR(3,4) = ZERO
  383.       CB2BAR(3,6) = ZERO
  384. C
  385.       CB2BAR(4,1) = ZERO
  386.       CB2BAR(4,2) = ZERO
  387.       CB2BAR(4,3) = ZERO
  388.       CB2BAR(4,4) = ZERO
  389. C
  390.       CB2BAR(5,4) = ZERO
  391.       CB2BAR(5,5) = ZERO
  392.       CB2BAR(5,6) = 1.D0
  393. C
  394.       CB2BAR(6,2) = ZERO
  395.       CB2BAR(6,5) = ZERO
  396.       CB2BAR(6,6) = 1.D0
  397. C
  398. C  ON A BESOIN DES VECTEURS GAMMA(8, ALPHA=1 A 4) = GS(8,4)
  399. C
  400.       DO I=1,24
  401.          XELOC(I) = XE(I)
  402.       END DO
  403. C
  404. C REMARQUE: ATTENTION, RR, MATRICE DU CORROTATIONNEL EST RANGEE PAR LIGNES!
  405. C
  406.       CALL CUBROT(XELOC,RR2)
  407.       CALL VECTRO(RR2,XELOC,1)
  408.       CALL CUBEBB(XELOC,B,VOL)
  409. C
  410. C CALCUL DES FONCTIONS DE FORME  GS
  411. C XXGB = X  * GB
  412. C
  413.       DO J=1,3
  414.          DO IA=1,4
  415.             XXGB(J,IA) = HOUXGB(XELOC(J),IA)
  416.          END DO
  417.       END DO
  418. C
  419. C GS = (BBB)  * XXGB
  420. C
  421.       DO I=1,8
  422.          DO J=1,4
  423.             GS(I,J)=0.D0
  424.          END DO
  425.       END DO
  426.       DO J=1,3
  427.          DO IA=1,4
  428.             DO I=1,8
  429.                GS(I,IA)=GS(I,IA) + B(J,I)*XXGB(J,IA)
  430.             END DO
  431.          END DO
  432.       END DO
  433. C
  434. C GS = GB - GS
  435. C
  436.       DO I=1,4
  437.          DO J=1,8
  438.             GS(J,I)= (GB(J,I) - GS(J,I))*UNS8
  439.          END DO
  440.       END DO
  441. C
  442. C CALCUL DE XXVB = X * VB
  443. C
  444.       XXVB(1)= - XELOC(1) + XELOC(4) + XELOC(7) - XELOC(10)
  445.      &      - XELOC(13) + XELOC(16) + XELOC(19) - XELOC(22)
  446.       XXVB(2)= - XELOC(2) - XELOC(5) + XELOC(8) + XELOC(11)
  447.      &      - XELOC(14) - XELOC(17) + XELOC(20) + XELOC(23)
  448.       XXVB(3)= - XELOC(3) - XELOC(6) - XELOC(9) - XELOC(12)
  449.      &      + XELOC(15) + XELOC(18) + XELOC(21) + XELOC(24)
  450. C
  451. C CALCUL DES RELATIONS CONTRAINTES ET DEFORMATIONS GENERALISEES
  452. C
  453.       HIJ(1) = UNS3*XXVB(2)*XXVB(3)/XXVB(1)
  454.       HIJ(2) = UNS3*XXVB(1)*XXVB(3)/XXVB(2)
  455.       HIJ(3) = UNS3*XXVB(2)*XXVB(1)/XXVB(3)
  456.       HIJ(4) = UNS3*XXVB(3)
  457.       HIJ(5) = UNS3*XXVB(1)
  458.       HIJ(6) = UNS3*XXVB(2)
  459. C
  460. C CALCUL DES COEFS A METTRE DANS KIJ POUR COMPOSER KSTAB
  461. C
  462. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  463. C ANCIENNE PROG A FARID (PLEXUS):
  464. C      XK11_1=(XLAMBDA+2*XMU)*HIJ(1)+XMU*HIJ(2)
  465. C      XK11_2=UNS3*((XLAMBDA+2*XMU)*HIJ(1)+XMU*(HIJ(2)+HIJ(3)))
  466. C      XK22_1=(XLAMBDA+2*XMU)*HIJ(2)+XMU*HIJ(1)
  467. C      XK22_2=UNS3*((XLAMBDA+2*XMU)*HIJ(2)+XMU*(HIJ(1)+HIJ(3)))
  468. C      XK33_1=XMU*(HIJ(1)+HIJ(2))
  469. C      XK33_2=UNS3*(PROPEL(1)*HIJ(3)+XMU*(HIJ(1)+HIJ(2)))
  470. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  471. C
  472. C   ICI IL FAUT PRENDRE LA MOYENNE DES MODULES D'YOUNG TANGENTS
  473. C
  474. C POUR LE SHB8 PLASTIQUE ET NON_LINEAIRE, METTRE IETAN A 1 DANS RIGID.FF
  475. C SINON METTRE IETAN A 0
  476. C
  477.       IF (IETAN.EQ.1) THEN
  478.          XYOUNGT = (EYG(1)+EYG(2)+EYG(3)+EYG(4)+EYG(5))/5.
  479.       ELSE
  480.          XYOUNGT = PROPEL(1)
  481.       ENDIF
  482. C
  483. C POUR RESOUDRE LE CISAILLEMENT TRANSVERSE:
  484. C
  485.       XYOUNGT = XYOUNGT*0.5*(COELA1+COELA2)
  486.       XLAMBDA = XYOUNGT*PROPEL(2)/(1-PROPEL(2)*PROPEL(2))
  487.       XMU     = 0.5 * XYOUNGT / (1+PROPEL(2))
  488.       XLAMBDA2= XLAMBDA + 2.D0*XMU
  489. C  ON MET LES LIGNES SUIVANTES EN COMMENTAIRE
  490. C
  491. C      XK11_1 = (XLAMBDA+2.D0*XMU)*HIJ(1)
  492. C      XK11_2 = UNS3*((XLAMBDA+2.D0*XMU)*HIJ(1))
  493. C      XK22_1 = (XLAMBDA+2.D0*XMU)*HIJ(2)
  494. C      XK22_2 = UNS3*((XLAMBDA+2.D0*XMU)*HIJ(2))
  495. C      XK33_1 = 0.
  496. C
  497.       XK11_1 = (CB2BAR(1,1)*(XLAMBDA2*CB2BAR(1,1) + XLAMBDA*CB2BAR(2,1))
  498.      &+ CB2BAR(2,1)*(XLAMBDA*CB2BAR(1,1) + XLAMBDA2*CB2BAR(2,1))
  499.      &+ XYOUNGT*CB2BAR(3,1)**2)*HIJ(1) + XMU*HIJ(2)*CB2BAR(4,1)**2
  500.       XK11_2 = (CB2BAR(1,2)*(XLAMBDA2*CB2BAR(1,2) + XLAMBDA*CB2BAR(2,2))
  501.      &+ CB2BAR(2,2)*(XLAMBDA*CB2BAR(1,2) + XLAMBDA2*CB2BAR(2,2))
  502. C     &+ XYOUNGT*CB2BAR(4,2)**2)*HIJ(1)*UNS3
  503.      &+ XYOUNGT*CB2BAR(3,2)**2)*HIJ(1)*UNS3
  504.      &+ XMU*UNS3*(HIJ(2)*CB2BAR(4,2)**2 + HIJ(3)*CB2BAR(6,2)**2)
  505. C
  506.       XK22_1 = (CB2BAR(1,3)*(XLAMBDA2*CB2BAR(1,3) + XLAMBDA*CB2BAR(2,3))
  507.      &+ CB2BAR(2,3)*(XLAMBDA*CB2BAR(1,3) + XLAMBDA2*CB2BAR(2,3))
  508.      &+ XYOUNGT*CB2BAR(3,3)**2)*HIJ(2) + XMU*HIJ(1)*CB2BAR(4,3)**2
  509.       XK22_2 = (CB2BAR(1,4)*(XLAMBDA2*CB2BAR(1,4) + XLAMBDA*CB2BAR(2,4))
  510.      &+ CB2BAR(2,4)*(XLAMBDA*CB2BAR(1,4) + XLAMBDA2*CB2BAR(2,4))
  511.      &+ XYOUNGT*CB2BAR(3,4)**2)*HIJ(2)*UNS3
  512.      &+ XMU*UNS3*(HIJ(1)*CB2BAR(4,4)**2 + HIJ(3)*CB2BAR(5,4)**2)
  513. C
  514. C      XK33_1 = XMU*(HIJ(2)*CB2BAR(5,5)**2 + HIJ(1)*CB2BAR(6,5))
  515.       XK33_1 = XMU*(HIJ(2)*CB2BAR(5,5)**2 + HIJ(1)*CB2BAR(6,5)**2)
  516.       XK33_2 = (CB2BAR(1,6)*(XLAMBDA2*CB2BAR(1,6) + XLAMBDA*CB2BAR(2,6))
  517.      &+ CB2BAR(2,6)*(XLAMBDA*CB2BAR(1,6) + XLAMBDA2*CB2BAR(2,6))
  518.      &+ XYOUNGT*CB2BAR(3,6)**2)*UNS3*HIJ(3)
  519.      &+ XMU*UNS3*(HIJ(2)*CB2BAR(5,6)**2 + HIJ(1)*CB2BAR(6,6)**2)
  520. C
  521. C ON MET LES SUIVANTES EN COMMENTAIRE POUR ANNULER LES TERMES CROISES
  522. C
  523. C      XK13_1 = XYOUNGT*CB2BAR(3,1)*HIJ(6)
  524. C      XK13_2 = XMU*CB2BAR(6,5)*HIJ(6)
  525. C      XK23_1 = XYOUNGT*CB2BAR(3,3)*HIJ(5)
  526. C      XK23_2 = XMU*CB2BAR(5,5)*HIJ(5)
  527. C      XK31_1 = XMU*CB2BAR(6,5)*HIJ(6)
  528. C      XK31_2 = XYOUNGT*CB2BAR(3,1)*HIJ(6)
  529. C      XK32_1 = XMU*CB2BAR(5,5)*HIJ(5)
  530. C      XK32_2 = XYOUNGT*CB2BAR(3,3)*HIJ(5)
  531. C
  532.       XK13_1 = ZERO
  533.       XK13_2 = ZERO
  534.       XK23_1 = ZERO
  535.       XK23_2 = ZERO
  536.       XK31_1 = ZERO
  537.       XK31_2 = ZERO
  538.       XK32_1 = ZERO
  539.       XK32_2 = ZERO
  540. C
  541. C  ON MET LES LIGNES SUIVANTES EN COMMENTAIRE
  542. C
  543. C      IF (IRDC.EQ.1) THEN
  544. C COMPORTEMENT SHB8 PLEXUS
  545. C         XK33_2 = XYOUNGT*HIJ(3)*UNS3
  546. C      ENDIF
  547. C
  548.       IF (IRDC.EQ.2) THEN
  549. C COMPORTEMENT C.P.
  550.          XK33_2 = 0.
  551.       ENDIF
  552. C
  553.       IF (IRDC.EQ.3) THEN
  554. C COMPORTEMENT LOI TRIDIM MMC 3D
  555.          XK11_1 = XCOOEF*(1-XNU)*HIJ(1)
  556.          XK11_2 = UNS3*(XCOOEF*(1-XNU)*HIJ(1))
  557.          XK22_1 = XCOOEF*(1-XNU)*HIJ(2)
  558.          XK22_2 = UNS3*(XCOOEF*(1-XNU)*HIJ(2))
  559.          XK33_1 = 0.
  560.          XK33_2 = XCOOEF*(1-XNU)*HIJ(3)*UNS3
  561.       ENDIF
  562. C
  563.       CALL ZDANUL(XK11,64)
  564.       CALL ZDANUL(XK22,64)
  565.       CALL ZDANUL(XK12,64)
  566.       CALL ZDANUL(XK21,64)
  567.       CALL ZDANUL(XK13,64)
  568.       CALL ZDANUL(XK23,64)
  569.       CALL ZDANUL(XK31,64)
  570.       CALL ZDANUL(XK32,64)
  571.       CALL ZDANUL(XK33,64)
  572. C
  573.       DO J=1,8
  574.          DO I=1,8
  575. C
  576. C IL FAUT CALCULER K11 K22 K33 K13 K23 K31 K32 MATRICES 8*8
  577. C
  578.       XK11(I,J)=XK11_1*GS(I,3)*GS(J,3)+XK11_2*GS(I,4)*GS(J,4)
  579.       XK22(I,J)=XK22_1*GS(I,3)*GS(J,3)+XK22_2*GS(I,4)*GS(J,4)
  580.       XK33(I,J)=XK33_1*GS(I,3)*GS(J,3)+XK33_2*GS(I,4)*GS(J,4)
  581. C
  582.       XK13(I,J)=XK13_1*GS(I,3)*GS(J,1)+XK13_2*GS(I,1)*GS(J,3)
  583.       XK23(I,J)=XK23_1*GS(I,3)*GS(J,2)+XK23_2*GS(I,2)*GS(J,3)
  584.       XK31(I,J)=XK31_1*GS(I,3)*GS(J,1)+XK31_2*GS(I,1)*GS(J,3)
  585.       XK32(I,J)=XK32_1*GS(I,3)*GS(J,2)+XK32_2*GS(I,2)*GS(J,3)
  586.  
  587.          END DO
  588.       END DO
  589. C
  590. C ASSEMBLAGE DE KSTAB
  591. C
  592.       CALL ASKSTB(XKSTAB,XK11,XK22,XK33,XK12,XK21,XK13,XK23,
  593.      &                XK31,XK32)
  594. C
  595. C
  596. C REMISE EN ORDRE DE KSTAB
  597. C
  598.       DO J=1,8
  599.          DO I=1,24
  600.             TMPKE(I,(J-1)*3+1)=XKSTAB(I,J)
  601.             TMPKE(I,(J-1)*3+2)=XKSTAB(I,J+8)
  602.             TMPKE(I,(J-1)*3+3)=XKSTAB(I,J+16)
  603.          END DO
  604.       END DO
  605.       CALL ZDANUL(XKSTAB,576)
  606.       DO I=1,8
  607.          DO J=1,24
  608.             XKSTAB((I-1)*3+1,J)=TMPKE(I,J)
  609.             XKSTAB((I-1)*3+2,J)=TMPKE(I+8,J)
  610.             XKSTAB((I-1)*3+3,J)=TMPKE(I+16,J)
  611.          END DO
  612.       END DO
  613. C
  614. C IL FAUT REPASSER DANS LE REPERE GLOBAL AVEC RR2^T . KSTAB . RR2
  615. C EN FAIT C'EST RR2. KSTAB .RR2^T CAR RR2 RANGEE PAR LIGNES
  616. C
  617.       CALL ASKSGL(XKSTAB,RR2)
  618. C
  619. C AJOUT DE KSTAB A KE
  620. C
  621.       DO J=1,24
  622.          DO I=1,24
  623.             RE(I,J)=RE(I,J)+XKSTAB(I,J)
  624.          END DO
  625.       END DO
  626.       ENDIF
  627. C
  628.       IF(ICLE.EQ.7.or.icle.eq.11)THEN
  629. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  630. C    icle=11 calcul de epsi                                                               C
  631. C ON CALCULE LES CONTRAINTES : SORTIE DANS OUT(30)                 C
  632. C                       CONTRAINTES LOCALES DANS CHAQUE COUCHE     C
  633. C                       SUR LA CONFIGURATION 1                     C
  634. C  LE DEPLACEMENT NODAL A L'AIR D'ETRE DANS WORK(1 A 24)           C
  635. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  636. C
  637. C UE: INCREMENT DE DEPLACEMENT NODAL, REPERE GLOBAL
  638. C
  639. C XE: DEBUT DU PAS
  640.       irdc=3
  641.       DO J=1,8
  642.          DO I=1,3
  643.             UE(I,J)=WORK((J-1)*3+I)
  644.          END DO
  645.       END DO
  646.       IGRDDEP = nint(D(1))
  647.       LAG     = nint(PROPEL(3))
  648. C
  649. C ON DEFINIT CMATLOC LOI MODIFIEE SHB8
  650. C
  651.       if(icle.eq.7) then
  652.       CALL ZDANUL(CMATLOC,36)
  653.       XLAMBDA       = PROPEL(1)*PROPEL(2)/(1-PROPEL(2)*PROPEL(2))
  654.       XMU             = 0.5 * PROPEL(1) / (1+PROPEL(2))
  655.       CMATLOC(1,1) = XLAMBDA+2*XMU
  656.       CMATLOC(2,2) = XLAMBDA+2*XMU
  657.       IF (IRDC.EQ.1) THEN
  658. C COMPORTEMENT SHB8 PLEXUS
  659.          CMATLOC(3,3) = PROPEL(1)
  660.       ENDIF
  661. C
  662.       IF (IRDC.EQ.2) THEN
  663. C COMPORTEMENT C.P.
  664.          CMATLOC(3,3) = 0
  665.       ENDIF
  666. C
  667.       CMATLOC(1,2) = XLAMBDA
  668.       CMATLOC(2,1) = XLAMBDA
  669.       CMATLOC(4,4) = XMU
  670.       CMATLOC(5,5) = XMU
  671.       CMATLOC(6,6) = XMU
  672. C
  673.       IF (IRDC.EQ.3) THEN
  674. C COMPORTEMENT LOI TRIDIM MMC 3D
  675.          XNU          = PROPEL(2)
  676.          XCOOEF       = PROPEL(1)/((1+XNU)*(1-2*XNU))
  677.          CMATLOC(1,1) = (1-XNU)*XCOOEF
  678.          CMATLOC(2,2) = (1-XNU)*XCOOEF
  679.          CMATLOC(3,3) = (1-XNU)*XCOOEF
  680.          CMATLOC(1,2) = XNU*XCOOEF
  681.          CMATLOC(2,1) = XNU*XCOOEF
  682.          CMATLOC(1,3) = XNU*XCOOEF
  683.          CMATLOC(3,1) = XNU*XCOOEF
  684.          CMATLOC(2,3) = XNU*XCOOEF
  685.          CMATLOC(3,2) = XNU*XCOOEF
  686.          CMATLOC(4,4) = (1-2*XNU)*0.5*XCOOEF
  687.          CMATLOC(5,5) = (1-2*XNU)*0.5*XCOOEF
  688.          CMATLOC(6,6) = (1-2*XNU)*0.5*XCOOEF
  689.       ENDIF
  690.       endif
  691. C
  692. C CALCUL DE BKSIP(3,8,IP) DANS REPERE DE REFERENCE
  693. C      BKSIP(1,*,IP) = VECTEUR BX AU POINT GAUSS IP
  694. C      BKSIP(2,*,IP) = VECTEUR BY AU POINT GAUSS IP
  695. C      BKSIP(3,*,IP) = VECTEUR BZ AU POINT GAUSS IP
  696. C
  697.       CALL CUBKSI(5,XXG5,BKSIP)
  698. C
  699.       DO IP=1,5
  700. C
  701. C DEFINITION DES 4 POINTS  COQUES
  702. C
  703.          ZETA  = XXG5(IP)
  704.          ZLAMB = 0.5*(1.-ZETA)
  705.          DO J=1,3
  706.             DO I=1,4
  707.                XCOQ(J,I) = ZLAMB*XE((I-1)*3+J)
  708.      &                    + (1.-ZLAMB)*XE((I-1+4)*3+J)
  709.             END DO
  710.          END DO
  711. C
  712. C CALCUL DE PPP 3*3 PASSAGE DE GLOBAL A LOCAL,COQUE
  713. C XCENT : COORD GLOBAL DU CENTRE DE L'ELEMENT
  714. C
  715.          CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,XBID)
  716. C
  717. C CALCUL DE B : U_GLOBAL ---> EPS_GLOBAL
  718. C
  719.          CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC)
  720. C
  721. C CALCUL DE EPS DANS LE REPERE GLOBAL: 1 POUR DEFORMATIONS LINEAIRES
  722. C                                     2 POUR TERMES CARRES EN PLUS
  723.          DO I=1,6
  724.             DEPS(I)=0.
  725.          ENDDO
  726.          IF (LAG.EQ.1) THEN
  727. C ON AJOUTE LA PARTIE NON-LINEAIRE DE EPS
  728.             CALL DSDX3D(2,B,UE,DEPS,DUSDX,8)
  729.          ELSE
  730.             CALL DSDX3D(1,B,UE,DEPS,DUSDX,8)
  731.          ENDIF
  732. C
  733. C SORTIE DE DUSDX DANS PROPEL(1 A 9 * 5 PT DE GAUSS)
  734. C POUR UTILISATION ULTERIEURE DANS Q8PKCN_SHB8
  735. C
  736.          CALL AEQBT(PPPT,PPP,3,3)
  737.          RR12(1,1) = DUSDX(1)
  738.          RR12(1,2) = DUSDX(2)
  739.          RR12(1,3) = DUSDX(3)
  740.          RR12(2,1) = DUSDX(4)
  741.          RR12(2,2) = DUSDX(5)
  742.          RR12(2,3) = DUSDX(6)
  743.          RR12(3,1) = DUSDX(7)
  744.          RR12(3,2) = DUSDX(8)
  745.          RR12(3,3) = DUSDX(9)
  746. *         CALL MULMAT(3,3,3,PPPT,RR12,RR2)
  747.          CALL MULMAT(RR2,PPPT,RR12,3,3,3)
  748. *         CALL MULMAT(3,3,3,RR2,PPP,RR12)
  749.          CALL MULMAT(RR12,RR2,PPP,3,3,3)
  750.          DUSDX(1)  = RR12(1,1)
  751.          DUSDX(2)  = RR12(1,2)
  752.          DUSDX(3)  = RR12(1,3)
  753.          DUSDX(4)  = RR12(2,1)
  754.          DUSDX(5)  = RR12(2,2)
  755.          DUSDX(6)  = RR12(2,3)
  756.          DUSDX(7)  = RR12(3,1)
  757.          DUSDX(8)  = RR12(3,2)
  758.          DUSDX(9)  = RR12(3,3)
  759.          DO I=1,9
  760.             PROPEL(I+(IP-1)*9)=DUSDX(I)
  761.          ENDDO
  762.          DO I=1,6
  763.             DEPS_LOC(I) = 0.
  764.             SIG_LOC(I)  = 0.
  765.          ENDDO
  766.          CALL CHRP3D(PPP,DEPS,DEPS_LOC,2)
  767. C
  768. C CALCUL DE SIGMA DANS LE REPERE LOCAL
  769. C
  770. *         CALL MULMAT(6,6,1,CMATLOC,DEPS_LOC,SIG_LOC)
  771.          if(icle.eq.7) then
  772.             CALL MULMAT(SIG_LOC,CMATLOC,DEPS_LOC,6,1,6)
  773.          endif
  774. C
  775. C CONTRAINTES ECRITES SOUS LA FORME:
  776. C               [SIG] = [S_11, S_22, S_33, S_12, S_23, S_13]
  777.          if(icle.eq.7) then
  778.          DO I=1,6
  779. C ON LAISSE LES CONTRAINTES DANS LE REPERE LOCAL POUR LA PLASTICITE
  780.             OUT((IP-1)*6+I)=SIG_LOC(I)
  781.          END DO
  782.          elseif(icle.eq.11)then
  783.          DO I=1,6
  784. C ON LAISSE LES CONTRAINTES DANS LE REPERE LOCAL POUR LA PLASTICITE
  785.             OUT((IP-1)*6+I)=DEPS_LOC(I)
  786.          END DO
  787.          endif
  788.       END DO
  789.       ENDIF
  790.       IF(icle.eq.12) then
  791. C propel(3)=1 on veut epsi à partir de sigma
  792. C propel(3)=2  on veut sigma à partir de epsi
  793. C en entrée work contient soit epsilon soit sigma
  794. C  work(1,2 ... contrainte 1 du point 1, contrainte 2 sdu point 1 ..
  795. C propel(1)=young
  796. C propel(2)=nu
  797. C propel(3) =1 si on veut des contraintes
  798. C           =2 si on veut des epsilons
  799. C
  800. C out en sortie out(1,2 .. epsi d1 du point 1 , epsi 2 du point 1 ..
  801.       irdc=3
  802.       icas=nint(propel(3))
  803. C
  804. C ON DEFINIT CMATLOC LOI MODIFIEE SHB8
  805. C
  806.       CALL ZDANUL(CMATLOC,36)
  807.       XLAMBDA       = PROPEL(1)*PROPEL(2)/(1-PROPEL(2)*PROPEL(2))
  808.       XMU             = 0.5 * PROPEL(1) / (1+PROPEL(2))
  809.       CMATLOC(1,1) = XLAMBDA+2*XMU
  810.       CMATLOC(2,2) = XLAMBDA+2*XMU
  811.       IF (IRDC.EQ.1) THEN
  812. C COMPORTEMENT SHB8 PLEXUS
  813.          CMATLOC(3,3) = PROPEL(1)
  814.       ENDIF
  815. C
  816.       IF (IRDC.EQ.2) THEN
  817. C COMPORTEMENT C.P.
  818.          CMATLOC(3,3) = 0
  819.       ENDIF
  820. C
  821.       CMATLOC(1,2) = XLAMBDA
  822.       CMATLOC(2,1) = XLAMBDA
  823.       CMATLOC(4,4) = XMU
  824.       CMATLOC(5,5) = XMU
  825.       CMATLOC(6,6) = XMU
  826. C
  827.       IF (IRDC.EQ.3) THEN
  828. C COMPORTEMENT LOI TRIDIM MMC 3D
  829.          XNU          = PROPEL(2)
  830.          XCOOEF       = PROPEL(1)/((1+XNU)*(1-2*XNU))
  831.          CMATLOC(1,1) = (1-XNU)*XCOOEF
  832.          CMATLOC(2,2) = (1-XNU)*XCOOEF
  833.          CMATLOC(3,3) = (1-XNU)*XCOOEF
  834.          CMATLOC(1,2) = XNU*XCOOEF
  835.          CMATLOC(2,1) = XNU*XCOOEF
  836.          CMATLOC(1,3) = XNU*XCOOEF
  837.          CMATLOC(3,1) = XNU*XCOOEF
  838.          CMATLOC(2,3) = XNU*XCOOEF
  839.          CMATLOC(3,2) = XNU*XCOOEF
  840.          CMATLOC(4,4) = (1-2*XNU)*0.5*XCOOEF
  841.          CMATLOC(5,5) = (1-2*XNU)*0.5*XCOOEF
  842.          CMATLOC(6,6) = (1-2*XNU)*0.5*XCOOEF
  843.       ENDIF
  844.       if(icas.eq.2) then
  845.         call invalm(cmatloc,6,3 , ier,1.d-12)
  846.         CMATLOC(4,4) =1.d0/CMATLOC(4,4)
  847.         CMATLOC(5,5) =1.d0/CMATLOC(5,5)
  848.         CMATLOC(6,6) =1.d0/CMATLOC(6,6)
  849.       endif
  850.       write(6,*) ' cmatloc '
  851.       write(6,*) (cmatloc(1,iop),iop=1,6)
  852.       write(6,*) (cmatloc(2,iop),iop=1,6)
  853.       write(6,*) (cmatloc(3,iop),iop=1,6)
  854.       write(6,*) (cmatloc(4,iop),iop=1,6)
  855.       write(6,*) (cmatloc(5,iop),iop=1,6)
  856.       write(6,*) (cmatloc(6,iop),iop=1,6)
  857.       write(6,*) 'epsi '
  858. C
  859.       DO IP=1,5
  860. C
  861. C DEFINITION DES 4 POINTS  COQUES
  862. C
  863.        write(6,*) (work((ip-1)*6 +iop),iop=1,6)
  864. C
  865. C CALCUL DE SIGMA DANS LE REPERE LOCAL
  866. C
  867.            CALL MULMAT(out((ip-1)*6+1),CMATLOC,work((ip-1)*6+1),6,1,6)
  868.       END DO
  869.       ENDIF
  870.  
  871.       IF(ICLE.EQ.8)THEN
  872. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  873. C                                                                  C
  874. C ON CALCULE BSIGMA: SORTIE DANS OUT(24)                           C
  875. C                    ENTREE DE SIGMA DANS WORK(DIM=30)             C
  876. C                    ENTREE DU MATERIAU DANS D(1) ET D(2)          C
  877. C                    ENTREE DE QIALPHA PAS PRECEDENT               C
  878. C                                      DANS RE(1 A 12)             C
  879. C                              QIALPHA(J,I)=RE(J+(I-1)*3)          C
  880. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  881.       LAG = nint(D(3))
  882.       IRDC=3
  883.       CALL ZDANUL(XKSIGTMP2,64)
  884.       DO J=1,8
  885.          DO I=1,3
  886.             F(I,J) = 0.D0
  887.          END DO
  888.       END DO
  889. C
  890. C CALCUL DE BKSIP(3,8,IP) DANS REPERE DE REFERENCE
  891. C      BKSIP(1,*,IP) = VECTEUR BX AU POINT GAUSS IP
  892. C      BKSIP(2,*,IP) = VECTEUR BY AU POINT GAUSS IP
  893. C      BKSIP(3,*,IP) = VECTEUR BZ AU POINT GAUSS IP
  894. C
  895.       CALL CUBKSI(5,XXG5,BKSIP)
  896.       DO IP=1,5
  897. C
  898. C RECHERCHE DE SIGMA DU POINT DE GAUSS GLOBAL
  899. C
  900.          DO I=1,6
  901. C
  902. C C'EST LES CONTRAINTES LOCALES POUR POUVOIR TRAITER LA PLASTICITE AVANT
  903. C LES CONTRAINTES SONT PASSEES DANS LA CONFI. A LA FIN DU PAS DANS Q8PKCN
  904. C
  905.             SIG_LOC(I)=WORK((IP-1)*6+I)
  906.          END DO
  907.          ZETA  = XXG5(IP)
  908.          ZLAMB = 0.5*(1.-ZETA)
  909.          DO J=1,3
  910.             DO I=1,4
  911.                XCOQ(J,I) = ZLAMB*XE((I-1)*3+J)
  912.      &                   + (1.-ZLAMB)*XE((I-1+4)*3+J)
  913.             END DO
  914.          END DO
  915.          CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,xBID)
  916. C
  917. C PASSAGE DES CONTRAINTES AU REPERE GLOBAL
  918. C
  919.          CALL CHRP3D(PPP,SIG_LOC,SIGMAG,1)
  920. C
  921.          CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC)
  922. C
  923. C CALCUL DE BQ.SIGMA SI LAGRANGIEN TOTAL
  924. C
  925.          IF (LAG.EQ.1) THEN
  926.              CALL ASKSIG(SIGMAG,B,XKSIGTMP1)
  927.              DO J=1,8
  928.                 DO I=1,8
  929.                    XKSIGTMP2(I,J) = XKSIGTMP2(I,J)
  930.      &                        + 4D0*AJAC*PXG5(IP)*XKSIGTMP1(I,J)
  931.                 ENDDO
  932.              ENDDO
  933.          ENDIF
  934. C
  935. C CALCUL DE B.SIGMA EN GLOBAL
  936. C
  937.          POIDS = 4D0 * AJAC * PXG5(IP)
  938.              DO K=1,8
  939.                 F(1,K) = F(1,K) + POIDS *
  940.      &          (B(1,K)*SIGMAG(1)+B(2,K)*SIGMAG(4)+B(3,K)*SIGMAG(6))
  941.                 F(2,K) = F(2,K) + POIDS *
  942.      &          (B(1,K)*SIGMAG(4)+B(2,K)*SIGMAG(2)+B(3,K)*SIGMAG(5))
  943.                 F(3,K) = F(3,K) + POIDS  *
  944.      &          (B(1,K)*SIGMAG(6)+B(2,K)*SIGMAG(5)+B(3,K)*SIGMAG(3))
  945.              END DO
  946. C
  947. C SI LAGRANGIEN TOTAL: AJOUT DE FQ A F
  948. C
  949.       END DO
  950.       IF (LAG.EQ.1) THEN
  951.          CALL ZDANUL(TMPKE,576)
  952.          DO KK=1,3
  953.             DO I=1,8
  954.                DO J=1,8
  955.                   TMPKE(I+(KK-1)*8,J+(KK-1)*8) = XKSIGTMP2(I,J)
  956.                ENDDO
  957.             ENDDO
  958.          ENDDO
  959.          CALL ZDANUL(TMPKE2,576)
  960.          DO J=1,8
  961.             DO I=1,24
  962.                TMPKE2(I,(J-1)*3+1)=TMPKE(I,J)
  963.                TMPKE2(I,(J-1)*3+2)=TMPKE(I,J+8)
  964.                TMPKE2(I,(J-1)*3+3)=TMPKE(I,J+16)
  965.             END DO
  966.          END DO
  967.          CALL ZDANUL(TMPKE,576)
  968.          DO I=1,8
  969.             DO J=1,24
  970.                TMPKE((I-1)*3+1,J)=TMPKE2(I,J)
  971.                TMPKE((I-1)*3+2,J)=TMPKE2(I+8,J)
  972.                TMPKE((I-1)*3+3,J)=TMPKE2(I+16,J)
  973.             END DO
  974.          END DO
  975. *         CALL MULMAT(24,24,1,TMPKE,PROPEL,FQ)
  976.          CALL MULMAT(FQ,TMPKE,PROPEL,24,1,24)
  977.          DO K=1,8
  978.             F(1,K) = F(1,K) + FQ((K-1)*3+1)
  979.             F(2,K) = F(2,K) + FQ((K-1)*3+2)
  980.             F(3,K) = F(3,K) + FQ((K-1)*3+3)
  981.          END DO
  982.       ENDIF
  983. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  984. C
  985. C ON CALCULE LA STABILISATION POUR LE CALCUL DE B.SIGMA
  986. C UE EST DANS PROPEL!
  987. C ON A BESOIN DU MATERIAU AUSSI!
  988. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  989. C
  990.       XYOUNG  = D(1)
  991.       XNU     = D(2)
  992.       XLAMBDA = XYOUNG*XNU/(1-XNU*XNU)
  993.       XMU     = 0.5 * XYOUNG / (1+XNU)
  994. C
  995. C DEPLACEMENT NODAL, REPERE GLOBAL : DANS PROPEL(1 A 24)
  996. C
  997.       DO I=1,24
  998.          XE_LOC(I)   = XE(I)
  999.          XE_LOCP(I)  = XE(I)-PROPEL(I)
  1000.          XE_LOC12(I) = XE(I)-0.5*PROPEL(I)
  1001.       ENDDO
  1002. C
  1003. C REMARQUE: ATTENTION, RR, MATRICE DU CORROTATIONNEL EST RANGEE PAR LIGNES!
  1004. C
  1005.       CALL CUBROT(XE_LOC,RR2)
  1006.       CALL CUBROT(XE_LOCP,RR1)
  1007.       CALL CUBROT(XE_LOC12,RR12)
  1008.       CALL VECTRO(RR2,XE_LOC,1)
  1009.       CALL VECTRO(RR1,XE_LOCP,1)
  1010.       CALL VECTRO(RR12,XE_LOC12,1)
  1011. C
  1012. C  ON A BESOIN DES VECTEURS GAMMA(8, ALPHA=1 A 4) = GS(8,4)
  1013. C  ORIGINE DE XE_LOC12 PAS BONNE, MAIS PAS GRAVE : DIFFERENTIELLES DANS B
  1014. C
  1015.       CALL CUBEBB(XE_LOC12,B,VOL)
  1016. C
  1017. C
  1018. C CALCUL DEPLACEMENT DEFORMANT DANS REPERE 1/2 PAS:
  1019. C
  1020.       DO I=1,24
  1021.          UDEF(I)=XE_LOC(I)-XE_LOCP(I)
  1022.       ENDDO
  1023. C
  1024. C CALCUL DES FONCTIONS DE FORME  GS
  1025. C XXGB = X  * GB
  1026. C
  1027.       DO IA=1,4
  1028.          DO J=1,3
  1029.             XXGB(J,IA) = HOUXGB(XE_LOC12(J),IA)
  1030.          END DO
  1031.       END DO
  1032. C
  1033. C GS = (BBB)  * XXGB
  1034. C
  1035.       DO J=1,4
  1036.          DO I=1,8
  1037.             GS(I,J)=0.D0
  1038.          END DO
  1039.       END DO
  1040.       DO J=1,3
  1041.          DO IA=1,4
  1042.             DO I=1,8
  1043.                GS(I,IA)=GS(I,IA) + B(J,I)*XXGB(J,IA)
  1044.             END DO
  1045.          END DO
  1046.       END DO
  1047. C
  1048. C GS = GB - GS
  1049. C
  1050.       DO I=1,4
  1051.          DO J=1,8
  1052.             GS(J,I)= (GB(J,I) - GS(J,I))*UNS8
  1053.          END DO
  1054.       END DO
  1055. C
  1056. C CALCUL DE XXVB = X * VB
  1057. C
  1058.       XXVB(1)= - XE_LOC12(1) + XE_LOC12(4) + XE_LOC12(7) - XE_LOC12(10)
  1059.      &      - XE_LOC12(13) + XE_LOC12(16) + XE_LOC12(19) - XE_LOC12(22)
  1060.       XXVB(2)= - XE_LOC12(2) - XE_LOC12(5) + XE_LOC12(8) + XE_LOC12(11)
  1061.      &      - XE_LOC12(14) - XE_LOC12(17) + XE_LOC12(20) + XE_LOC12(23)
  1062.       XXVB(3)= - XE_LOC12(3) - XE_LOC12(6) - XE_LOC12(9) - XE_LOC12(12)
  1063.      &      + XE_LOC12(15) + XE_LOC12(18) + XE_LOC12(21) + XE_LOC12(24)
  1064. C
  1065. C CALCUL DES RELATIONS CONTRAINTES ET DEFORMATIONS GENERALISEES
  1066. C
  1067.       HIJ(1) = UNS3*XXVB(2)*XXVB(3) / XXVB(1)
  1068.       HIJ(2) = UNS3*XXVB(1)*XXVB(3) / XXVB(2)
  1069.       HIJ(3) = UNS3*XXVB(2)*XXVB(1) / XXVB(3)
  1070.       HIJ(4) = UNS3*XXVB(3)
  1071.       HIJ(5) = UNS3*XXVB(1)
  1072.       HIJ(6) = UNS3*XXVB(2)
  1073. C
  1074. C CALCUL DES DEFORMATIONS GENERALISEES
  1075. C
  1076.       DO IA=1,4
  1077.          DO J=1,3
  1078.             AUX=0.D0
  1079.             DO I=1,8
  1080.                AUX = AUX + GS(I,IA)*UDEF((I-1)*3+J)
  1081.             END DO
  1082.             PQIALPH(J,IA) = AUX
  1083.          END DO
  1084.       END DO
  1085. C
  1086. C CALCUL DES CONTRAINTES GENERALISEES
  1087. C
  1088.       DO I=1,4
  1089.          DO J=1,3
  1090.             QIALPHA(J,I)=RE(J+(I-1)*3,1)
  1091.          END DO
  1092.       END DO
  1093. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1094. C              PROGRAMMATION FARID (PLEXUS)                        C
  1095. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1096. C       QIALPHA(1,1) = QIALPHA(1,1)
  1097. C       QIALPHA(2,2) = QIALPHA(2,2)
  1098. C       QIALPHA(3,3) = QIALPHA(3,3) + XMU*((HIJ(2)+HIJ(1))*PQIALPH(3,3))
  1099. C       QIALPHA(1,2) = QIALPHA(1,2)
  1100. C       QIALPHA(2,3) = QIALPHA(2,3) + ((XLAMBDA+2.D0*XMU)*HIJ(2)
  1101. C      &           + XMU*HIJ(1))*PQIALPH(2,3)
  1102. C       QIALPHA(1,3) = QIALPHA(1,3) + ((XLAMBDA+2.D0*XMU)*HIJ(1)
  1103. C      &           + XMU*HIJ(2))*PQIALPH(1,3)
  1104. C       QIALPHA(2,1) = QIALPHA(2,1)
  1105. C       QIALPHA(3,2) = QIALPHA(3,2)
  1106. C       QIALPHA(3,1) = QIALPHA(3,1)
  1107. C       QIALPHA(1,4) = QIALPHA(1,4) + UNS3*((XLAMBDA+2.D0*XMU)*HIJ(1)
  1108. C      &          + XMU*(HIJ(2)+HIJ(3)))*PQIALPH(1,4)
  1109. C       QIALPHA(2,4) = QIALPHA(2,4) + UNS3*((XLAMBDA+2.D0*XMU)*HIJ(2)
  1110. C      &           + XMU*(HIJ(1)+HIJ(3)))*PQIALPH(2,4)
  1111. C       QIALPHA(3,4) = QIALPHA(3,4) + UNS3*(XYOUNG*HIJ(3)
  1112. C     &          + XMU*(HIJ(1)+HIJ(2)))*PQIALPH(3,4)
  1113. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1114.       IF (IRDC.EQ.1) THEN
  1115.          QIALPHA(1,1) = QIALPHA(1,1)
  1116.          QIALPHA(2,2) = QIALPHA(2,2)
  1117.          QIALPHA(3,3) = QIALPHA(3,3)
  1118.          QIALPHA(1,2) = QIALPHA(1,2)
  1119.          QIALPHA(2,3) = QIALPHA(2,3)
  1120.      &                + ((XLAMBDA+2.D0*XMU)*HIJ(2))*PQIALPH(2,3)
  1121.          QIALPHA(1,3) = QIALPHA(1,3)
  1122.      &                + ((XLAMBDA+2.D0*XMU)*HIJ(1))*PQIALPH(1,3)
  1123.          QIALPHA(2,1) = QIALPHA(2,1)
  1124.          QIALPHA(3,2) = QIALPHA(3,2)
  1125.          QIALPHA(3,1) = QIALPHA(3,1)
  1126.          QIALPHA(1,4) = QIALPHA(1,4)
  1127.      &                + UNS3*((XLAMBDA+2*XMU)*HIJ(1))*PQIALPH(1,4)
  1128.          QIALPHA(2,4) = QIALPHA(2,4)
  1129.      &                + UNS3*((XLAMBDA+2*XMU)*HIJ(2))*PQIALPH(2,4)
  1130. C PROJECTION LEGAY POUR LE TERME QIALPHA(3,4)
  1131. C COMPORTEMENT SHB8 PLEXUS
  1132. CC         QIALPHA(3,4) = QIALPHA(3,4) + XYOUNG*HIJ(3)*UNS3*PQIALPH(3,4)
  1133. C PROJECTION 28 LE TERME QIALPHA(3,4) EST MODIFIE
  1134.          QIALPHA(3,4)=QIALPHA(3,4)+UNS3*(HIJ(1)+HIJ(2))*XMU*PQIALPH(3,4)
  1135.       ENDIF
  1136. C A MODIFIER POUR LA PROJECTION 28 (CAS DES C.P.)
  1137.       IF (IRDC.EQ.2) THEN
  1138.          QIALPHA(1,1) = QIALPHA(1,1)
  1139.          QIALPHA(2,2) = QIALPHA(2,2)
  1140.          QIALPHA(3,3) = QIALPHA(3,3)
  1141.          QIALPHA(1,2) = QIALPHA(1,2)
  1142.          QIALPHA(2,3) = QIALPHA(2,3)
  1143.      &                + ((XLAMBDA+2.D0*XMU)*HIJ(2))*PQIALPH(2,3)
  1144.          QIALPHA(1,3) = QIALPHA(1,3)
  1145.      &                + ((XLAMBDA+2.D0*XMU)*HIJ(1))*PQIALPH(1,3)
  1146.          QIALPHA(2,1) = QIALPHA(2,1)
  1147.          QIALPHA(3,2) = QIALPHA(3,2)
  1148.          QIALPHA(3,1) = QIALPHA(3,1)
  1149.          QIALPHA(1,4) = QIALPHA(1,4)
  1150.      &                + UNS3*((XLAMBDA+2*XMU)*HIJ(1))*PQIALPH(1,4)
  1151.          QIALPHA(2,4) = QIALPHA(2,4)
  1152.      &                + UNS3*((XLAMBDA+2*XMU)*HIJ(2))*PQIALPH(2,4)
  1153. C COMPORTEMENT C.P.
  1154.          QIALPHA(3,4) = QIALPHA(3,4)
  1155.       ENDIF
  1156.       IF (IRDC.EQ.3) THEN
  1157. C COMPORTEMENT LOI TRIDIM MMC 3D
  1158. C A MODIFIER POUR LA PROJECTION 28 (CAS DE LA LOI 3D)
  1159.          XCOOEF       = XYOUNG/((1+XNU)*(1-2*XNU))
  1160.          QIALPHA(1,1) = QIALPHA(1,1)
  1161.          QIALPHA(2,2) = QIALPHA(2,2)
  1162.          QIALPHA(3,3) = QIALPHA(3,3)
  1163.          QIALPHA(1,2) = QIALPHA(1,2)
  1164.          QIALPHA(2,3) = QIALPHA(2,3)
  1165.      &                + ((XCOOEF*(1-XNU))*HIJ(2))*PQIALPH(2,3)
  1166.          QIALPHA(1,3) = QIALPHA(1,3)
  1167.      &                + ((XCOOEF*(1-XNU))*HIJ(1))*PQIALPH(1,3)
  1168.          QIALPHA(2,1) = QIALPHA(2,1)
  1169.          QIALPHA(3,2) = QIALPHA(3,2)
  1170.          QIALPHA(3,1) = QIALPHA(3,1)
  1171.          QIALPHA(1,4) = QIALPHA(1,4)
  1172.      &                + UNS3*((XCOOEF*(1-XNU))*HIJ(1))*PQIALPH(1,4)
  1173.          QIALPHA(2,4) = QIALPHA(2,4)
  1174.      &                + UNS3*((XCOOEF*(1-XNU))*HIJ(2))*PQIALPH(2,4)
  1175.          QIALPHA(3,4) = QIALPHA(3,4)
  1176.      &                + XCOOEF*(1-XNU)*HIJ(3)*UNS3*PQIALPH(3,4)
  1177.       ENDIF
  1178. CCCCCCCCCCCCCCCCCCC
  1179. C
  1180. C SAUVEGARDE DES FORCES DE STABILISATION
  1181. C
  1182.       DO I=1,4
  1183.          DO J=1,3
  1184.             RE(J+(I-1)*3,1)=QIALPHA(J,I)
  1185.          END DO
  1186.       END DO
  1187. C
  1188. C CALCUL DES FORCES DE HOURGLASS FHG
  1189. C
  1190.       DO I=1,8
  1191.          DO J=1,3
  1192.             FHG(J,I)=0.D0
  1193.          ENDDO
  1194.       END DO
  1195.       DO J=1,3
  1196.          DO I=1,8
  1197.             DO IA=1,4
  1198.                FHG(J,I) = FHG(J,I) + QIALPHA(J,IA)*GS(I,IA)
  1199.             END DO
  1200.          END DO
  1201.       END DO
  1202. C ON REPASSE AU REPERE GLOBAL
  1203.       DO I=1,3
  1204.          DO J=1,8
  1205.             FHG24((J-1)*3+I)=FHG(I,J)
  1206.          ENDDO
  1207.       ENDDO
  1208.       CALL VECTRO(RR12,FHG24,2)
  1209. C
  1210. C AJOUT DE LA STABILISATION AU B SIGMA DEJA CALCULE
  1211. C
  1212.       DO J=1,3
  1213.          DO I=1,8
  1214.            F(J,I) = F(J,I) + FHG24((I-1)*3+J)
  1215.          ENDDO
  1216.       ENDDO
  1217. C
  1218. C ATTENTION A L'ORDRE DE OUT
  1219. C
  1220.       DO I=1,3
  1221.          DO J=1,8
  1222.             OUT((J-1)*3+I) = F(I,J)
  1223.          END DO
  1224.       END DO
  1225.       ENDIF
  1226.       IF(ICLE.EQ.5)THEN
  1227. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1228. C                                                                  C
  1229. C ON CALCULE LES FORCES DE PRESSION                                C
  1230. C                                                                  C
  1231. C      ENTREES: PROPEL(1)   : VALEUR DE LA PRESSION                C
  1232. C                 XE          : GEOMETRIE ELEMENT                  C
  1233. C           WORK(1 A 3) : DIRECTION DE PRESSION                    C
  1234. C               PROPEL(2)   : 0 RECHERCHE DE LA FACE               C
  1235. C                             1 OU 2 : FACE INF. OU SUP. (POUR     C
  1236. C                                       PRESSION SUIVEUSE)         C
  1237. C     SORTIE: OUT FORCE DE PRESSION                                C
  1238. C               WORK(4)     : NUMERO DE FACE POUR LA PRESSION      C
  1239. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1240.       DO I=1,8
  1241.          DO J=1,3
  1242.             F(J,I)=0.
  1243.          END DO
  1244.       END DO
  1245. C
  1246. C LECTURE: PRESSION A METTRE AUX 4 NOEUDS, FACE ET DIRECTION
  1247. C
  1248.       PRESS=PROPEL(1)
  1249.       IFACE=nint(PROPEL(2))
  1250.       IF(IFACE.EQ.0)THEN
  1251.          DO I=1,3
  1252.             DIREC(I)=WORK(I)
  1253.          END DO
  1254.       ENDIF
  1255. C
  1256. C COORD. DES POINTS MILIEUX
  1257. C
  1258.       DO J=1,3
  1259.          X56(J) =( XE((5-1)*3+J)+XE((6-1)*3+J) )/2.
  1260.          X67(J) =( XE((6-1)*3+J)+XE((7-1)*3+J) )/2.
  1261.          X78(J) =( XE((7-1)*3+J)+XE((8-1)*3+J) )/2.
  1262.          X85(J) =( XE((8-1)*3+J)+XE((5-1)*3+J) )/2.
  1263.          X12(J) =( XE((1-1)*3+J)+XE((2-1)*3+J) )/2.
  1264.          X23(J) =( XE((2-1)*3+J)+XE((3-1)*3+J) )/2.
  1265.          X34(J) =( XE((3-1)*3+J)+XE((4-1)*3+J) )/2.
  1266.          X41(J) =( XE((4-1)*3+J)+XE((1-1)*3+J) )/2.
  1267.       END DO
  1268. C
  1269. C NORMALE RENTRANTE FACE SUP: V56_78 ^ V85_67
  1270. C
  1271.       IF(IFACE.EQ.0.OR.IFACE.EQ.2)THEN
  1272.          DO J=1,3
  1273.             V56_78(J)=X78(J)-X56(J)
  1274.             V85_67(J)=X67(J)-X85(J)
  1275.          END DO
  1276.          VSUP(1) = V56_78(2) * V85_67(3) - V56_78(3) * V85_67(2)
  1277.          VSUP(2) = V56_78(3) * V85_67(1) - V56_78(1) * V85_67(3)
  1278.          VSUP(3) = V56_78(1) * V85_67(2) - V56_78(2) * V85_67(1)
  1279.       ENDIF
  1280. C
  1281. C NORMALE RENTRANTE FACE INF: V41_23 ^ V12_34
  1282. C
  1283.       IF(IFACE.EQ.0.OR.IFACE.EQ.1)THEN
  1284.          DO J=1,3
  1285.             V12_34(J)=X34(J)-X12(J)
  1286.             V41_23(J)=X23(J)-X41(J)
  1287.          END DO
  1288.          VINF(1) = V41_23(2) * V12_34(3) - V41_23(3) * V12_34(2)
  1289.          VINF(2) = V41_23(3) * V12_34(1) - V41_23(1) * V12_34(3)
  1290.          VINF(3) = V41_23(1) * V12_34(2) - V41_23(2) * V12_34(1)
  1291.       ENDIF
  1292.       IF(IFACE.EQ.0) THEN
  1293. C
  1294. C RECHERCHE DE LA BONNE FACE: PREMIER PASSAGE
  1295. C
  1296.          SCALSUP = 0.
  1297.          SCALINF = 0.
  1298.          DO I=1,3
  1299.             SCALSUP = SCALSUP+VSUP(I)*DIREC(I)
  1300.             SCALINF = SCALINF+VINF(I)*DIREC(I)
  1301.          END DO
  1302. C
  1303. C TEST ERREUR
  1304. C
  1305.          IF((SCALSUP*SCALINF).GE.0)THEN
  1306.             IF ((SCALSUP*SCALINF).EQ.0)THEN
  1307.                WRITE(6,*)'**      ERREUR SUR LES SURFACES        **'
  1308.                WRITE(6,*)'**       DE SHB8 POUR PRESSION         **'
  1309.                WRITE(6,*)'**    PRODUIT SCALAIRE ÉGAL A ZERO     **'
  1310.                WRITE(6,*)'**  PROBLÈME DE DIRECTION DE PRESSION  **'
  1311.             ENDIF
  1312.             IF ((SCALSUP*SCALINF).GT.0)THEN
  1313.                WRITE(6,*)'**    ERREUR SUR LES SURFACES   **'
  1314.                WRITE(6,*)'**     DE SHB8 POUR PRESSION    **'
  1315.                WRITE(6,*)'**  LES 2 FACES SONT POSSIBLES  **'
  1316.             ENDIF
  1317.             STOP
  1318.          ENDIF
  1319. C
  1320. C SI SCAL <0 : PRESSION SUR CETTE FACE
  1321. C
  1322.          IF(SCALSUP.LT.0)THEN
  1323.             WORK(4) = 1
  1324.             IFACE   = 1
  1325.          ENDIF
  1326.          IF(SCALINF.LT.0)THEN
  1327.             WORK(4) = 2
  1328.             IFACE   = 2
  1329.          ENDIF
  1330.       ENDIF
  1331.       IF(IFACE.EQ.2)THEN
  1332. C
  1333. C BOUCLE SUR LES 4 NOEUDS DE LA FACE SUPERIEURE
  1334. C
  1335.          DO I=5,8
  1336.             DO J=1,3
  1337.                F(J,I)=  PRESS * VSUP(J) / 4.
  1338.             END DO
  1339.          END DO
  1340.       ENDIF
  1341.       IF(IFACE.EQ.1)THEN
  1342. C
  1343. C BOUCLE SUR LES 4 NOEUDS DE LA FACE INFERIEURE
  1344. C
  1345.          DO I=1,4
  1346.             DO J=1,3
  1347.                F(J,I)=  PRESS * VINF(J) / 4.
  1348.             END DO
  1349.          END DO
  1350.       ENDIF
  1351. C
  1352. C ATTENTION A L'ORDRE DE OUT
  1353. C
  1354.       DO J=1,8
  1355.          DO I=1,3
  1356.             OUT((J-1)*3+I)=F(I,J)
  1357.          END DO
  1358.       END DO
  1359.       ENDIF
  1360.       IF(ICLE.EQ.9)THEN
  1361. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1362. C CALCUL DE K.SIGMA                                                C
  1363. C          EN ENTREE DANS WORK : SIGMA NORMALEMENT LONGUEUR 30     C
  1364. C                         PROPEL(1): 1 POUR RE=RE+KSIGMA           C
  1365. C                         PROPEL(1): 0 POUR RE=KSIGMA              C
  1366. C             SORTIE           : RE(24*24)=RE(24*24)+KSIGMA        C
  1367. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1368. C
  1369. C CALCUL DE B (1 2 3) AUX 5 POINTS DE GAUSS
  1370. C
  1371.       CALL CUBKSI(5,XXG5,BKSIP)
  1372.       DO J=1,8
  1373.          DO I=1,8
  1374.           XKSIGTMP2(I,J)=0.
  1375.          ENDDO
  1376.       ENDDO
  1377. C
  1378. C
  1379. C CALCUL DE LA CONTRAINTE MOYENNE SUR LES 5 POINTS DE GAUSS POUR LE
  1380. C KSIGMA DE STABILISATION CALCULE PLUS LOIN
  1381. C
  1382.       DO J=1,6
  1383.           SIGMOY(J)=0.D0
  1384.       ENDDO
  1385. C
  1386.       DO J=1,6
  1387.          DO IP=1,5
  1388.            SIGMOY(J)=SIGMOY(J)+WORK((IP-1)*6+J)
  1389.          ENDDO
  1390.            SIGMOY(J)=UNS5*SIGMOY(J)
  1391.       ENDDO
  1392. C
  1393. C DEBUT DE LA BOUCLE SUR LES 5 PTS GAUSS
  1394. C
  1395.       DO IP=1,5
  1396. C
  1397. C CALCUL DE B
  1398. C
  1399.          CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC)
  1400. C
  1401. C CALCUL DE MATRICE DE PASSAGE POUR POUVOIR CALCULER LES CONTRAINTES DANS
  1402. C LE REPERE GLOBAL
  1403. C
  1404.          DO I=1,6
  1405. C C'EST LES CONTRAINTES LOCALES POUR POUVOIR TRAITER LA PLASTICITE AVANT
  1406.             SIG_LOC(I)=WORK((IP-1)*6+I)
  1407.          END DO
  1408.          ZETA  = XXG5(IP)
  1409.          ZLAMB = 0.5*(1.-ZETA)
  1410.          DO J=1,3
  1411.             DO I=1,4
  1412.                XCOQ(J,I) = ZLAMB*XE((I-1)*3+J)
  1413.      &                     + (1.-ZLAMB)*XE((I-1+4)*3+J)
  1414.             END DO
  1415.          END DO
  1416.          CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,XBID)
  1417. C
  1418. C PASSAGE DES CONTRAINTES AU REPERE GLOBAL
  1419. C
  1420.          CALL CHRP3D(PPP,SIG_LOC,SIGMAG,1)
  1421.          CALL ASKSIG(SIGMAG,B,XKSIGTMP1)
  1422.          DO J=1,8
  1423.             DO I=1,8
  1424.                XKSIGTMP2(I,J) = XKSIGTMP2(I,J)
  1425.      &                          + 4D0*AJAC*PXG5(IP)*XKSIGTMP1(I,J)
  1426.             ENDDO
  1427.          ENDDO
  1428.       ENDDO
  1429.       CALL ZDANUL(TMPKE,576)
  1430.       DO KK=1,3
  1431.          DO I=1,8
  1432.             DO J=1,8
  1433.                TMPKE(I+(KK-1)*8,J+(KK-1)*8) = XKSIGTMP2(I,J)
  1434.             ENDDO
  1435.          ENDDO
  1436.       ENDDO
  1437. C
  1438. C ON MET DE L'ORDRE:
  1439. C
  1440.       CALL ZDANUL(TMPKE2,576)
  1441.       DO J=1,8
  1442.          DO I=1,24
  1443.             TMPKE2(I,(J-1)*3+1)=TMPKE(I,J)
  1444.             TMPKE2(I,(J-1)*3+2)=TMPKE(I,J+8)
  1445.             TMPKE2(I,(J-1)*3+3)=TMPKE(I,J+16)
  1446.          END DO
  1447.       END DO
  1448.       CALL ZDANUL(TMPKE,576)
  1449.       DO I=1,8
  1450.          DO J=1,24
  1451.             TMPKE((I-1)*3+1,J)=TMPKE2(I,J)
  1452.             TMPKE((I-1)*3+2,J)=TMPKE2(I+8,J)
  1453.             TMPKE((I-1)*3+3,J)=TMPKE2(I+16,J)
  1454.          END DO
  1455.       END DO
  1456.       IPROPEL=nint(PROPEL(1))
  1457.       IF(IPROPEL.EQ.0) CALL SHIFTD(TMPKE,RE,576)
  1458.       IF(IPROPEL.EQ.1) CALL AEQAPB(RE,TMPKE,576)
  1459. C
  1460. C
  1461. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1462. C                                                                              C
  1463. C MATRICE DE STABILISATION POUR K.SIGMA: PAS DE BOUCLE SUR LES POINTS DE GAUSS C
  1464. C                                                                              C
  1465. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1466. C
  1467. C
  1468. C
  1469. C  ON A BESOIN DES VECTEURS GAMMA(8, ALPHA=1 A 4) = GS(8,4)
  1470. C
  1471. CCC      DO I=1,24
  1472. CCC         XELOC(I) = XE(I)
  1473. CCC      END DO
  1474. C
  1475. C REMARQUE: ATTENTION, RR, MATRICE DU CORROTATIONNEL EST RANGEE PAR LIGNES!
  1476. C
  1477. CCC      CALL CUBEROT(XELOC,RR2)
  1478. CCC      CALL VECTROT(RR2,XELOC,1)
  1479. CCC      CALL CUBE_BB(XELOC,B,VOL)
  1480. C
  1481. C CALCUL DES FONCTIONS DE FORME  GS
  1482. C XXGB = X  * GB
  1483. C
  1484. CCC      DO J=1,3
  1485. CCC         DO IA=1,4
  1486. CCC            XXGB(J,IA) = HOUXGB(XELOC(J),IA)
  1487. CCC         END DO
  1488. CCC      END DO
  1489. C
  1490. C GS = (BBB)  * XXGB
  1491. C
  1492. CCC      DO I=1,8
  1493. CCC         DO J=1,4
  1494. CCC            GS(I,J)=0.D0
  1495. CCC         END DO
  1496. CCC      END DO
  1497. CCC      DO J=1,3
  1498. CCC         DO IA=1,4
  1499. CCC            DO I=1,8
  1500. CCC               GS(I,IA)=GS(I,IA) + B(J,I)*XXGB(J,IA)
  1501. CCC            END DO
  1502. CCC         END DO
  1503. CCC      END DO
  1504. C
  1505. C GS = GB - GS
  1506. C
  1507. CCC      DO I=1,4
  1508. CCC         DO J=1,8
  1509. CCC            GS(J,I)= (GB(J,I) - GS(J,I))*UNS8
  1510. CCC         END DO
  1511. CCC      END DO
  1512. C
  1513. C CALCUL DE XXVB = X * VB
  1514. C
  1515. CCC      XXVB(1)= - XELOC(1) + XELOC(4) + XELOC(7) - XELOC(10)
  1516. CCC     &      - XELOC(13) + XELOC(16) + XELOC(19) - XELOC(22)
  1517. CCC      XXVB(2)= - XELOC(2) - XELOC(5) + XELOC(8) + XELOC(11)
  1518. CCC     &      - XELOC(14) - XELOC(17) + XELOC(20) + XELOC(23)
  1519. CCC      XXVB(3)= - XELOC(3) - XELOC(6) - XELOC(9) - XELOC(12)
  1520. CCC     &      + XELOC(15) + XELOC(18) + XELOC(21) + XELOC(24)
  1521. C
  1522. C CALCUL DES RELATIONS CONTRAINTES ET DEFORMATIONS GENERALISEES
  1523. C
  1524. CCC      HIJ(1) = UNS3*XXVB(2)*XXVB(3)/XXVB(1)
  1525. CCC      HIJ(2) = UNS3*XXVB(1)*XXVB(3)/XXVB(2)
  1526. CCC      HIJ(3) = UNS3*XXVB(2)*XXVB(1)/XXVB(3)
  1527. CCC      HIJ(4) = UNS3*XXVB(3)
  1528. CCC      HIJ(5) = UNS3*XXVB(1)
  1529. CCC      HIJ(6) = UNS3*XXVB(2)
  1530. C
  1531. C CALCUL DES COEFS A METTRE DANS KIJ POUR COMPOSER KSTAB
  1532. C
  1533. C
  1534. C ICI IL FAUT PRENDRE LA MOYENNE DES CONTRAINTES SUR LES POINTS DE GAUSS
  1535. C CECI EST DEJA CALCULE PLUS HAUT DANS LA VARIABLE SIGMOY(6)
  1536. C
  1537. C
  1538. C  POUR LA PROJECTION 28 ON OBTIENT
  1539. C
  1540. CCC      XK11_1 = SIGMOY(1)*HIJ(1)
  1541. CCC      XK11_2 = SIGMOY(1)*HIJ(1)*UNS3
  1542. CCC      XK11_3 = SIGMOY(6)*HIJ(6)
  1543. CCC      XK11_4 = SIGMOY(6)*HIJ(6)
  1544. C
  1545. CCC      XK22_1 = SIGMOY(2)*HIJ(2)
  1546. CCC      XK22_2 = SIGMOY(2)*HIJ(2)*UNS3
  1547. CCC      XK22_3 = SIGMOY(5)*HIJ(5)
  1548. CCC      XK22_4 = SIGMOY(5)*HIJ(5)
  1549. C
  1550. CCC      XK33_1 = ZERO
  1551. CCC      XK33_2 = (SIGMOY(1)*HIJ(1)+SIGMOY(2)*HIJ(2))*UNS3
  1552. C
  1553. CCC      XK13_1 = ZERO
  1554. CCC      XK13_2 = ZERO
  1555. CCC      XK23_1 = ZERO
  1556. CCC      XK23_2 = ZERO
  1557. CCC      XK31_1 = ZERO
  1558. CCC      XK31_2 = ZERO
  1559. CCC      XK32_1 = ZERO
  1560. CCC      XK32_2 = ZERO
  1561. C
  1562. C
  1563. CCC      CALL ZDANUL(XK11,64)
  1564. CCC      CALL ZDANUL(XK22,64)
  1565. CCC      CALL ZDANUL(XK12,64)
  1566. CCC      CALL ZDANUL(XK21,64)
  1567. CCC      CALL ZDANUL(XK13,64)
  1568. CCC      CALL ZDANUL(XK23,64)
  1569. CCC      CALL ZDANUL(XK31,64)
  1570. CCC      CALL ZDANUL(XK32,64)
  1571. CCC      CALL ZDANUL(XK33,64)
  1572. C
  1573. CCC      DO J=1,8
  1574. CCC         DO I=1,8
  1575. C
  1576. C IL FAUT CALCULER K11 K22 K33 K13 K23 K31 K32 MATRICES 8*8
  1577. C
  1578. CCC      XK11(I,J)=XK11_1*GS(I,3)*GS(J,3)+XK11_2*GS(I,4)*GS(J,4)
  1579. C     &          +XK11_3*GS(I,3)*GS(J,1)+XK11_4*GS(I,1)*GS(J,3)
  1580. CCC      XK22(I,J)=XK22_1*GS(I,3)*GS(J,3)+XK22_2*GS(I,4)*GS(J,4)
  1581. C     &          +XK22_3*GS(I,3)*GS(J,2)+XK22_4*GS(I,2)*GS(J,3)
  1582. CCC      XK33(I,J)=XK33_1*GS(I,3)*GS(J,3)+XK33_2*GS(I,4)*GS(J,4)
  1583. C
  1584. CCC      XK13(I,J)=XK13_1*GS(I,3)*GS(J,1)+XK13_2*GS(I,1)*GS(J,3)
  1585. CCC      XK23(I,J)=XK23_1*GS(I,3)*GS(J,2)+XK23_2*GS(I,2)*GS(J,3)
  1586. CCC      XK31(I,J)=XK31_1*GS(I,3)*GS(J,1)+XK31_2*GS(I,1)*GS(J,3)
  1587. CCC      XK32(I,J)=XK32_1*GS(I,3)*GS(J,2)+XK32_2*GS(I,2)*GS(J,3)
  1588. C
  1589. CCC         END DO
  1590. CCC      END DO
  1591. C
  1592. C ASSEMBLAGE DE KSTAB
  1593. C
  1594. CCC      CALL ASSE_KSTAB(XKSTAB,XK11,XK22,XK33,XK12,XK21,XK13,XK23,
  1595. CCC     &                XK31,XK32)
  1596. C
  1597. C
  1598. C REMISE EN ORDRE DE KSTAB
  1599. C
  1600. CCC      DO J=1,8
  1601. CCC         DO I=1,24
  1602. CCC            TMPKE(I,(J-1)*3+1)=XKSTAB(I,J)
  1603. CCC            TMPKE(I,(J-1)*3+2)=XKSTAB(I,J+8)
  1604. CCC            TMPKE(I,(J-1)*3+3)=XKSTAB(I,J+16)
  1605. CCC         END DO
  1606. CCC      END DO
  1607. CCC      CALL ZDANUL(XKSTAB,576)
  1608. CCC      DO I=1,8
  1609. CCC         DO J=1,24
  1610. CCC            XKSTAB((I-1)*3+1,J)=TMPKE(I,J)
  1611. CCC            XKSTAB((I-1)*3+2,J)=TMPKE(I+8,J)
  1612. CCC            XKSTAB((I-1)*3+3,J)=TMPKE(I+16,J)
  1613. CCC         END DO
  1614. CCC      END DO
  1615. C
  1616. C IL FAUT REPASSER DANS LE REPERE GLOBAL AVEC RR2^T . KSTAB . RR2
  1617. C EN FAIT C'EST RR2. KSTAB .RR2^T CAR RR2 RANGEE PAR LIGNES
  1618. C
  1619. CCC      CALL ASSE_KSTAB_GL(XKSTAB,RR2)
  1620. C
  1621. C AJOUT DE KSTAB A KE
  1622. C
  1623. CCC      DO J=1,24
  1624. CCC         DO I=1,24
  1625. CCC            RE(I,J)=RE(I,J)+XKSTAB(I,J)
  1626. CCC         END DO
  1627. CCC      END DO
  1628. C
  1629. C
  1630. C
  1631.       ENDIF
  1632. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1633.       IF(ICLE.EQ.3)THEN
  1634. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1635. C                                                                  C
  1636. C   CALCUL DE LA MATRICE DE MASSE                                  C
  1637. C           ENTREE: PROPEL(1) = MASSE VOLUMIQUE                    C
  1638. C                 XE(1 A 14)= COORDONNEES DES NOEUDS               C
  1639. C           SORTIE: RE(24,24) = MATRICE MASSE                      C
  1640. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1641.       XRO = PROPEL(1)
  1642.       CALL ZDANUL(RE,576)
  1643. C
  1644. C CALCUL DU PLAN MOYEN:
  1645. C
  1646.       DO J=1,3
  1647.          DO I=1,4
  1648.             XCOQ(J,I) = 0.5*XE((I-1)*3+J)+0.5*XE((I-1+4)*3+J)
  1649.          END DO
  1650.       END DO
  1651. C CALCUL DE LA SURFACE MOYENNE:
  1652.       DO J=1,3
  1653.          X12(J) = (XCOQ(J,1)+XCOQ(J,2))*0.5
  1654.          X34(J) = (XCOQ(J,3)+XCOQ(J,4))*0.5
  1655.          X23(J) = (XCOQ(J,2)+XCOQ(J,3))*0.5
  1656.          X41(J) = (XCOQ(J,4)+XCOQ(J,1))*0.5
  1657.       END DO
  1658.       DO J=1,3
  1659.          V12_34(J) = X34(J)-X12(J)
  1660.          V41_23(J) = X23(J)-X41(J)
  1661.       END DO
  1662. C
  1663. C CALCUL DE V41_23 ^ V12_34:
  1664. C
  1665.       VINF(1) = V41_23(2) * V12_34(3) - V41_23(3) * V12_34(2)
  1666.       VINF(2) = V41_23(3) * V12_34(1) - V41_23(1) * V12_34(3)
  1667.       VINF(3) = V41_23(1) * V12_34(2) - V41_23(2) * V12_34(1)
  1668. C
  1669. C CALCUL DU VOLUME DE L'ELEMENT:
  1670. C
  1671.       CALL CUBKSI(5,XXG5,BKSIP)
  1672.       CALL CUCALB(BKSIP(1,1,3),XE,B,AJAC)
  1673.       XVOLUME  = AJAC*8.
  1674.       CALL ZDANUL(XMASSE,64)
  1675.       XXXMASSE = 0.
  1676.       DO I=1,8
  1677.          XMASSE(I,I) = (8./27.)*AJAC*XRO
  1678.          XXXMASSE    = XXXMASSE+XMASSE(I,I)
  1679.       ENDDO
  1680.       DO I=1,8
  1681.          DO J=1,8
  1682.             IF(I.NE.J)THEN
  1683.                XMASSE(I,J) = (1./64.)*(2.+VB(I,1)*VB(J,1)*2./3.)
  1684.      &                        *(2.+VB(I,2)*VB(J,2)*2./3.)
  1685.      &                        *(2.+VB(I,3)*VB(J,3)*2./3.)
  1686.      &                        *AJAC*XRO
  1687.                XXXMASSE    = XXXMASSE+XMASSE(I,I)
  1688.             ENDIF
  1689.          ENDDO
  1690.       ENDDO
  1691.       CALL ZDANUL(TMPKE,576)
  1692.       DO K=1,3
  1693.          DO I=1,8
  1694.             DO J=1,8
  1695.                TMPKE(I+(K-1)*8,J+(K-1)*8) = XMASSE(I,J)
  1696.             ENDDO
  1697.          ENDDO
  1698.       ENDDO
  1699.       DO J=1,8
  1700.          DO I=1,24
  1701.             TMPKE2(I,(J-1)*3+1)=TMPKE(I,J)
  1702.             TMPKE2(I,(J-1)*3+2)=TMPKE(I,J+8)
  1703.             TMPKE2(I,(J-1)*3+3)=TMPKE(I,J+16)
  1704.          END DO
  1705.       END DO
  1706.       CALL ZDANUL(RE,576)
  1707.       DO I=1,8
  1708.          DO J=1,24
  1709.             RE((I-1)*3+1,J)=TMPKE2(I,J)
  1710.             RE((I-1)*3+2,J)=TMPKE2(I+8,J)
  1711.             RE((I-1)*3+3,J)=TMPKE2(I+16,J)
  1712.          END DO
  1713.       END DO
  1714.       ENDIF
  1715.       IF(ICLE.EQ.10)THEN
  1716. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1717. C                                                                  C
  1718. C   CALCUL DE LA MATRICE KP                                        C
  1719. C                       ENTREE: PRESSION DANS PROPEL(1)            C
  1720. C                             KTAFL DANS PROPEL(2)                 C
  1721. C                       RE=RE+KP                                   C
  1722. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1723.       PRESSI = PROPEL(1)
  1724.       IFACE  = nint(PROPEL(2))
  1725.       ITEST  = 0
  1726.       IF (IFACE.EQ.1) THEN
  1727. C
  1728. C DEFINITION DE LA NORMALE N0, FACE INFERIEURE
  1729. C
  1730.          DO I=1,4
  1731.             DO J=1,3
  1732.                XCOQ(J,I) = XE((I-1)*3+J)
  1733.             END DO
  1734.          END DO
  1735.          DO J=1,3
  1736.             XCOQ14(J) = (XCOQ(J,1)+XCOQ(J,4))/2.
  1737.             XCOQ23(J) = (XCOQ(J,2)+XCOQ(J,3))/2.
  1738.             XCOQ34(J) = (XCOQ(J,3)+XCOQ(J,4))/2.
  1739.             XCOQ12(J) = (XCOQ(J,1)+XCOQ(J,2))/2.
  1740.          ENDDO
  1741.          DO J=1,3
  1742.             XCOQ12_34(J) = XCOQ34(J)-XCOQ12(J)
  1743.             XCOQ14_23(J) = XCOQ23(J)-XCOQ14(J)
  1744.          ENDDO
  1745. C
  1746.          XNORMALE(1) = XCOQ14_23(2)*XCOQ12_34(3)
  1747.      &                   -XCOQ14_23(3)*XCOQ12_34(2)
  1748.          XNORMALE(2) = XCOQ14_23(3)*XCOQ12_34(1)
  1749.      &                   -XCOQ14_23(1)*XCOQ12_34(3)
  1750.          XNORMALE(3) = XCOQ14_23(1)*XCOQ12_34(2)
  1751.      &                   -XCOQ14_23(2)*XCOQ12_34(1)
  1752.          XNORMALE(1) = -XNORMALE(1)
  1753.          XNORMALE(2) = -XNORMALE(2)
  1754.          XNORMALE(3) = -XNORMALE(3)
  1755. C
  1756.          XXG5(1) = -1.
  1757.          XXG5(2) = 0.
  1758.          XXG5(3) = 0.
  1759.          XXG5(4) = 0.
  1760.          XXG5(5) = 0.
  1761.          ITEST   = 1
  1762.       ENDIF
  1763.       IF (IFACE.EQ.2) THEN
  1764. C
  1765. C DEFINITION DE LA NORMALE N0, FACE SUPERIEURE
  1766. C
  1767.          DO I=1,4
  1768.             DO J=1,3
  1769.                XCOQ(J,I) = XE(12+(I-1)*3+J)
  1770.             END DO
  1771.          END DO
  1772.          DO J=1,3
  1773.             XCOQ14(J)=(XCOQ(J,1)+XCOQ(J,4))/2.
  1774.             XCOQ23(J)=(XCOQ(J,2)+XCOQ(J,3))/2.
  1775.             XCOQ34(J)=(XCOQ(J,3)+XCOQ(J,4))/2.
  1776.             XCOQ12(J)=(XCOQ(J,1)+XCOQ(J,2))/2.
  1777.          ENDDO
  1778.          DO J=1,3
  1779.             XCOQ12_34(J)=XCOQ34(J)-XCOQ12(J)
  1780.             XCOQ14_23(J)=XCOQ23(J)-XCOQ14(J)
  1781.          ENDDO
  1782. C
  1783.          XNORMALE(1)=XCOQ14_23(2)*XCOQ12_34(3)
  1784.      &                 -XCOQ14_23(3)*XCOQ12_34(2)
  1785.          XNORMALE(2)=XCOQ14_23(3)*XCOQ12_34(1)
  1786.      &                 -XCOQ14_23(1)*XCOQ12_34(3)
  1787.          XNORMALE(3)=XCOQ14_23(1)*XCOQ12_34(2)
  1788.      &                 -XCOQ14_23(2)*XCOQ12_34(1)
  1789. C
  1790.          XXG5(1)=+1.
  1791.          XXG5(2)=0.
  1792.          XXG5(3)=0.
  1793.          XXG5(4)=0.
  1794.          XXG5(5)=0.
  1795.          ITEST=1
  1796.       ENDIF
  1797.       IF(ITEST.EQ.0)THEN
  1798.          WRITE(6,*)'ARRET DANS SHB8, CALCUL DE KP, FACE NON DEFINIE'
  1799.          STOP
  1800.       ENDIF
  1801. C
  1802. C CALCUL DE B SUR LA FACE INF OU SUP, SUIVANT XXG5
  1803. C
  1804. C DANS UE: ON MET L'EQUIVALENT DE BKSIP, CALCULER
  1805. C                 POUR LES 4 NOEUDS DE LA FACE IFACE
  1806. C
  1807.       CALL CUKSKP(IFACE,UE)
  1808. C
  1809.       XNORM = SQRT(XNORMALE(1)*XNORMALE(1)+XNORMALE(2)*XNORMALE(2)
  1810.      &                                     +XNORMALE(3)*XNORMALE(3))
  1811.       DO J=1,3
  1812.          XNORMALE(J) = XNORMALE(J)/XNORM
  1813.       ENDDO
  1814.       CALL CUCBKP(IFACE,UE,XE,BKP,AJAC)
  1815.       CALL ASSKP2(IFACE,BKP,XNORMALE,XKP)
  1816. C
  1817. C REMISE EN ORDRE DE KP
  1818. C
  1819.       DO I=1,24
  1820.          DO J=1,8
  1821.             TMPKE(I,(J-1)*3+1)=XKP(I,J)
  1822.             TMPKE(I,(J-1)*3+2)=XKP(I,J+8)
  1823.             TMPKE(I,(J-1)*3+3)=XKP(I,J+16)
  1824.          END DO
  1825.       END DO
  1826.       CALL ZDANUL(XKP,576)
  1827.       DO I=1,8
  1828.          DO J=1,24
  1829.             XKP((I-1)*3+1,J)=TMPKE(I,J)
  1830.             XKP((I-1)*3+2,J)=TMPKE(I+8,J)
  1831.             XKP((I-1)*3+3,J)=TMPKE(I+16,J)
  1832.          END DO
  1833.       END DO
  1834.       DO I=1,24
  1835.          DO J=1,24
  1836.             RE(I,J)=RE(I,J)+XKP(I,J)*XNORM*PRESSI/8.
  1837.          ENDDO
  1838.       ENDDO
  1839.       ENDIF
  1840.       RETURN
  1841.       END
  1842.  
  1843.  
  1844.  
  1845.  
  1846.  
  1847.  
  1848.  
  1849.  

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