1 : C PB301 SOURCE MAGN 10/05/18 21:16:31 6675 2 : SUBROUTINE PB301(XREF,X,PG,FN,GR,FM,GM,ND,NP,MP,NG,NPG,NOM2,ITYPI) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C************************************************************************ 6 : C 7 : C CALCULE LES FONCTIONS DE FORME D'UN : SEG3 Element quadratique Q2 8 : C : SEG3 Element macro 2xQ1 9 : C (pression Q1 dans les 2 cas) 10 : C Cf VNIMP pression continue 11 : C ITYPI=0 1 pts d'integration de type Gauss (sous integration) 12 : C ITYPI=1 pts d'integration de type Gauss (en nombre suffisant) 13 : C ITYPI=2 pts d'integration de type Gauss Lobatto i.e. sommets element 14 : C 15 : C REMARQUE : 16 : C En coordonnées cylindriques si le point milieux (2) est 17 : C pile poil au milieu, la matrice masse diagonale s'annule sur l'axe. 18 : C /1 19 : C | 2 pi x N1 dx = 0. 20 : C /0 21 : C Pour remédier à cela on décale un peu le point mileux sur l'élément 22 : C de référence, alfa=0.50000000000005 au lieu de 0.5 . 23 : C Voir aussi pb902.eso et pb702.eso. 24 : C Si on rapprochait les points milieux de l'axe (alfa = 0.45 par exemple) 25 : C la matrice masse diagonale serait négative sur l'axe, ce que l'on ne 26 : C veut pas non plus. En conséquence le choix qu'on a fait doit être 27 : C accompagné d'une orientation correcte de l'élément courant. 28 : C************************************************************************ 29 : 30 : DIMENSION X(NPG),XREF(ND,NP) 31 : DIMENSION FN(NP,NPG),GR(ND,NP,NPG),PG(NPG) 32 : REAL*8 FM(MP,NPG),GM(ND,MP,NPG) 33 : DIMENSION U(5),H(5) 34 : CHARACTER*4 NOM2 35 : -INC CCOPTIO 36 : 37 : XREF(1,1)=0.D0 38 : XREF(1,2)=0.5D0 39 : XREF(1,3)=1.D0 40 : 41 : IF(IFOMOD.EQ.0)THEN 42 : alfa=0.50000000000005D0 43 : ELSE 44 : alfa=0.5D0 45 : ENDIF 46 : 47 : h1=1.D0/alfa 48 : h2=-1.D0/alfa/(alfa - 1.D0) 49 : h3=1.D0/(alfa - 1.D0) 49 : CALL CALUHG(U,H,NG) 50 : A=0.D0 51 : B=1.D0 52 : IF(ITYPI.EQ.2)THEN 53 : X(1)=0.D0 54 : X(2)=alfa 55 : X(3)=1.D0 56 : DO 4 L=1,3 57 : PG(L)=1.D0/3.D0 58 : 4 CONTINUE 59 : ELSE 59 : CALL CALG1(A,B,NG,H,U,X,PG) 60 : ENDIF 61 : 62 : C write(6,*)'Xg=',(X(I),I=1,NPG) 63 : C write(6,*)'Pg=',(Pg(I),I=1,NPG) 64 : 65 : IF(NOM2.NE.'MCF1')THEN 66 : DO 1 L=1,NPG 67 : FN(1,L)=h1*(X(L)-1.D0)*(X(L)-alfa) 68 : FN(2,L)=h2*X(L)*(1.D0-X(L)) 69 : FN(3,L)=h3*X(L)*(alfa-X(L)) 70 : GR(1,1,L)=h1*(2.D0*X(L)- 1.D0 - alfa) 71 : GR(1,2,L)=h2*(1.D0-2.D0*X(L)) 72 : GR(1,3,L)=h3*(alfa-2.D0*X(L)) 73 : 1 CONTINUE 74 : ELSE 75 : C Cas MACRO ELEMENT 76 : DO 10 L=1,NPG 77 : IF(X(L).LE.0.5D0)THEN 78 : FN(1,L)=1.D0-2.D0*X(L) 79 : FN(2,L)=2.D0*X(L) 80 : FN(3,L)=0.D0 81 : GR(1,1,L)=-2.D0 82 : GR(1,2,L)=2.D0 83 : GR(1,3,L)=0.D0 84 : ELSE 85 : FN(1,L)=0.D0 86 : FN(2,L)=2.D0-2.D0*X(L) 87 : FN(3,L)=2.D0*X(L)-1.D0 88 : GR(1,1,L)=0.D0 89 : GR(1,2,L)=-2.D0 90 : GR(1,3,L)=2.D0 91 : ENDIF 92 : 10 CONTINUE 93 : ENDIF 94 : 95 : IF(MP.EQ.1)THEN 96 : DO 2 L=1,NPG 97 : FM(1,L)=1.D0 98 : GM(1,1,L)=0.D0 99 : 2 CONTINUE 100 : ELSEIF(MP.EQ.2)THEN 101 : 102 : FM(1,1)= 1.D0 103 : FM(1,2)= 0.D0 104 : 105 : FM(2,1)= 0.D0 106 : FM(2,2)= 1.D0 107 : 107 : CALL INITD(GM,(ND*MP*NPG),0.D0) 108 : 109 : ENDIF 110 : IF(NOM2.EQ.'P1P1'.OR.NOM2.EQ.'PFP1'.OR.NOM2.EQ.'MCF1')THEN 111 : DO 3 L=1,NPG 112 : C 113 : FM(1,L)=1.D0-X(L) 114 : FM(2,L)=X(L) 115 : GM(1,1,L)=-1.D0 116 : GM(1,2,L)=1.D0 117 : 3 CONTINUE 118 : 119 : ENDIF 120 : 121 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 122 : C WRITE(6,101) 123 : C DO 104 L=1,NPG 124 : C write(6,*)'FN npg=',L 125 : C WRITE(6,1002)(FN(I,L),I=1,NP) 126 : C write(6,*)'GR(1 npg=',L 127 : C WRITE(6,1002)(GR(1,I,L),I=1,NP) 128 : C write(6,*)'GR(2 npg=',L 129 : C WRITE(6,1002)(GR(2,I,L),I=1,NP) 130 : C write(6,*)'FM npg=',L 131 : C WRITE(6,1002)(FM(I,L),I=1,MP) 132 : C WRITE(6,1002)FM 133 : C WRITE(6,1002)GM 134 : C104 CONTINUE 135 : C WRITE(6,101) 136 : RETURN 137 : 1002 FORMAT(10(1X,1PD11.4)) 138 : 1001 FORMAT(20(1X,I5)) 139 : 101 FORMAT(1X,'... SUB PB301 ... FN,GR ',9(10H..........)/) 140 : C 141 : END 142 : 143 : 144 : 145 : 146 : 147 : 148 : 149 : 150 : 151 :
© Cast3M 2003 - All rights reserved.
Disclaimer