1 : C JACOBI SOURCE CHAT 05/01/13 00:48:25 5004 2 : 3 : C======================================================================= 4 : C= J A C O B I = 5 : C= ----------- = 6 : C= Fonction : = 7 : C= ---------- = 8 : C= Calcul du jacobien (DJac) et des derivees des fonctions de forme = 9 : C= par rapport aux coordonnees REELLES x,y,z (SHPTOT(2 a 4,*)) en un = 10 : C= point donne d'un element (en general, un point de Gauss) = 11 : C= = 12 : C= Parametres : (E)=Entree (S)=Sortie = 13 : C= ------------ = 14 : C= XEL(3,NBNO) (E) Coordonnees GLOBALES des noeuds de l element = 15 : C= SHP(6,NBNO) (E/S) Valeurs des fonctions de forme et de leurs = 16 : C= derivees au point considere (point de Gauss) = 17 : C= IDIM (E) Dimension du probleme traite (1 a 3) = 18 : C= Cas particulier de l'element JOI3 : IDIM = 86 = 19 : C= NBNO (E) Nombre de NOEUDS de l element fini = 20 : C= DJac (S) Jacobien calcule en ce point de Gauss = 21 : C= (egal a 0 si probleme en ce point) = 22 : C= = 23 : C= Remarque : = 24 : C= ---------- = 25 : C= Lors de l'entree dans le sous-programme, SHP(2 a 4,*) contient = 26 : C= les DERIVEES des fonctions de forme par rapport aux coordonnees = 27 : C= de REFERENCE Qsi,Eta,Dzeta. = 28 : C= En sortie du sous-programme, SHP(2 a 4,*) contient les DERIVEES = 29 : C= des fonctions de FORME par rapport aux coordonnees REELLES x,y,z. = 30 : C======================================================================= 31 : 32 : SUBROUTINE JACOBI (XEL,SHP,IDIM,NBNO,DJac) 33 : 34 : IMPLICIT INTEGER(I-N) 35 : IMPLICIT REAL*8 (A-H,O-Z) 36 : 37 : DIMENSION XEL(3,*),SHP(6,*) 38 : 39 : PARAMETER (XZero=0.,XUn=1.) 40 : 41 : dJInv=XZero 42 : 43 : C =========================== 44 : C 1 - Cas de la DIMENSION 1 45 : C =========================== 46 : IF (IDIM.EQ.1) THEN 47 : DJac=XZero 48 : DO i=1,NBNO 49 : DJac=DJac+SHP(2,i)*XEL(1,i) 50 : ENDDO 51 : IF (DJac.NE.XZero) dJInv=XUn/DJac 52 : DO i=1,NBNO 53 : SHP(2,i)=SHP(2,i)*dJInv 54 : ENDDO 55 : C =========================== 56 : C 2 - Cas de la DIMENSION 2 57 : C =========================== 58 : ELSE IF (IDIM.EQ.2) THEN 59 : dXdQsi=XZero 60 : dXdEta=XZero 61 : dYdQsi=XZero 62 : dYdEta=XZero 63 : DO i=1,NBNO 64 : dXdQsi=dXdQsi+SHP(2,i)*XEL(1,i) 65 : dXdEta=dXdEta+SHP(3,i)*XEL(1,i) 66 : dYdQsi=dYdQsi+SHP(2,i)*XEL(2,i) 67 : dYdEta=dYdEta+SHP(3,i)*XEL(2,i) 68 : ENDDO 69 : DJac=dXdQsi*dYdEta-dXdEta*dYdQsi 70 : IF (DJac.NE.XZero) dJInv=XUn/DJac 71 : dQsidX= dYdEta*dJInv 72 : dQsidY=-dXdEta*dJInv 73 : dEtadX=-dYdQsi*dJInv 74 : dEtaDY= dXdQsi*dJInv 75 : DO i=1,NBNO 76 : V1=SHP(2,i)*dQsidX+SHP(3,i)*dEtadX 77 : SHP(3,i)=SHP(2,i)*dQsidY+SHP(3,i)*dEtadY 78 : SHP(2,i)=V1 79 : ENDDO 80 : C =========================== 81 : C 3 - Cas de la DIMENSION 3 82 : C =========================== 83 : ELSE IF (IDIM.EQ.3) THEN 84 : D11=XZero 85 : D21=XZero 86 : D31=XZero 87 : D12=XZero 88 : D22=XZero 89 : D32=XZero 90 : D13=XZero 91 : D23=XZero 92 : D33=XZero 93 : DO i=1,NBNO 94 : D11=D11+SHP(2,i)*XEL(1,i) 95 : D21=D21+SHP(3,i)*XEL(1,i) 96 : D31=D31+SHP(4,i)*XEL(1,i) 97 : D12=D12+SHP(2,i)*XEL(2,i) 98 : D22=D22+SHP(3,i)*XEL(2,i) 99 : D32=D32+SHP(4,i)*XEL(2,i) 100 : D13=D13+SHP(2,i)*XEL(3,i) 101 : D23=D23+SHP(3,i)*XEL(3,i) 102 : D33=D33+SHP(4,i)*XEL(3,i) 103 : ENDDO 104 : DInv11=D22*D33-D23*D32 105 : DInv12=D32*D13-D12*D33 106 : DInv13=D12*D23-D22*D13 107 : DJac=D11*DInv11+D21*DInv12+D31*DInv13 108 : IF (DJac.NE.XZero) dJInv=XUn/DJac 109 : DInv11=DInv11*dJInv 110 : DInv12=DInv12*dJInv 111 : DInv13=DInv13*dJInv 112 : DInv21=(D23*D31-D21*D33)*dJInv 113 : DInv22=(D11*D33-D13*D31)*dJInv 114 : DInv23=(D21*D13-D11*D23)*dJInv 115 : DInv31=(D21*D32-D22*D31)*dJInv 116 : DInv32=(D12*D31-D11*D32)*dJInv 117 : DInv33=(D11*D22-D12*D21)*dJInv 118 : DO i=1,NBNO 119 : V1=DInv11*SHP(2,i)+DInv12*SHP(3,i)+DInv13*SHP(4,i) 120 : V2=DInv21*SHP(2,i)+DInv22*SHP(3,i)+DInv23*SHP(4,i) 121 : V3=DInv31*SHP(2,i)+DInv32*SHP(3,i)+DInv33*SHP(4,i) 122 : SHP(2,i)=V1 123 : SHP(3,i)=V2 124 : SHP(4,i)=V3 125 : ENDDO 126 : C ======================================= 127 : C 4 - Cas particulier de l'element JOI3 128 : C ======================================= 129 : ELSE IF (IDIM.EQ.86) THEN 130 : dXdQsi=XZero 131 : dYdQsi=XZero 132 : DO i=1,NBNO 133 : dXdQsi=dXdQsi+SHP(2,i)*XEL(1,i) 134 : dYdQsi=dYdQsi+SHP(2,i)*XEL(2,i) 135 : ENDDO 136 : DJac=SQRT(dXdQsi*dXdQsi+dYdQsi*dYdQsi) 137 : ENDIF 138 : 139 : RETURN 140 : END 141 : 142 : 143 :
© Cast3M 2003 - All rights reserved.
Disclaimer