1 : C XLAPL SOURCE CHAT 05/01/13 04:14:37 5004 2 : SUBROUTINE XLAPL(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI, 3 : & COEF,COG,COES,KITT,KJTT,IK1, 4 : & LE,NBEL,K0,XCOOR,AIMPL,IKOMP, 5 : & AF1,AF2,AF3, 6 : & AS1,AS2,AS3, 7 : & NINC,IHV,IARG,S2) 8 : C 9 : IMPLICIT INTEGER(I-N) 10 : IMPLICIT REAL*8 (A-H,O-Z) 11 : C************************************************************************ 12 : C 13 : C CALCUL MATRICE ELEMENTAIRE DU LAPLACIEN 14 : C 15 : C 16 : C************************************************************************ 17 : 18 : C DIMENSION FN(NP,NPG),HR(KES,NP,NPG),PGSQ(NPG),RPG(NPG) 19 : DIMENSION LE(NP,NBEL) 20 : DIMENSION XYZ(IDIM,NP),FN(NP,NPG),GR(IDIM,NP,NPG),PG(NPG) 21 : DIMENSION HR(IDIM,NP,NPG),PGSQ(NPG),RPG(NPG) 22 : DIMENSION AF1(NBEL,NP,NP),AF2(NBEL,NP,NP),AF3(NBEL,NP,NP) 23 : DIMENSION AS1(NBEL,NP,NP),AS2(NBEL,NP,NP),AS3(NBEL,NP,NP) 24 : DIMENSION CL(9),COEF(1),COES(IDIM,IDIM,*) 25 : DIMENSION AF(4,4,2,2),XCOOR(*) 26 : C 27 : DIMENSION COE(3,3),S2(IDIM,1),COG(2,1) 28 : 29 : C write(6,*)' XLAPL KITT,KJTT,IK1 =',KITT,KJTT,IK1 30 : NK=K0 31 : DO 108 KE=1,NBEL 32 : NK=NK+1 33 : DO 109 I=1,NP 34 : J=LE(I,KE) 35 : DO 109 N=1,IDIM 36 : XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N) 37 : 109 CONTINUE 38 : 38 : CALL CALJBC(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES, 39 : *IDIM,NP,NPG,IAXI,AIRE) 40 : C 41 : IAX=0 42 : C WRITE(6,*) 43 : C WRITE(6,*)' KITT ',KITT,' KJTT ',KJTT 44 : C if(iax.eq.0)go to 108 45 : IF(IHV.EQ.1)IAX=IAXI 46 : IF(KITT.EQ.4)GO TO 70 47 : IF(KITT.EQ.3)GO TO 80 48 : IF(KITT.EQ.1)GO TO 80 49 : IF(KJTT.EQ.4)GO TO 5 50 : IC=1+(1-IK1)*(NK-1) 51 : C WRITE(6,*)' IC ',IC ,COEF(IC) 52 : DO 1 I=1,NP 53 : DO 1 J=1,NP 54 : U=0.D0 55 : DO 2 L=1,NPG 56 : V=0.D0 57 : DO 3 N=1,IDIM 58 : V=V+HR(N,I,L)*HR(N,J,L) 59 : 3 CONTINUE 60 : U=U+V*PGSQ(L) 61 : 2 CONTINUE 62 : U=U*COEF(IC) 63 : AF1(KE,J,I)=AF1(KE,J,I)+U*AIMPL 64 : AS1(KE,J,I)=AS1(KE,J,I)+U*(AIMPL-1.D0) 65 : 1 CONTINUE 66 : 67 : IF(NINC.GE.2)THEN 68 : DO 1762 I=1,NP 69 : DO 1762 J=1,NP 70 : AF2(KE,I,J)=AF1(KE,I,J) 71 : AS2(KE,I,J)=AS1(KE,I,J) 72 : 1762 CONTINUE 73 : ENDIF 74 : IF(NINC.GE.3)THEN 75 : DO 1763 I=1,NP 76 : DO 1763 J=1,NP 77 : AF3(KE,I,J)=AF1(KE,I,J) 78 : AS3(KE,I,J)=AS1(KE,I,J) 79 : 1763 CONTINUE 80 : ENDIF 81 : 82 : IF(IHV.EQ.1)THEN 83 : IF(IAXI.EQ.1)THEN 84 : DO 41 I=1,NP 85 : DO 41 J=1,NP 86 : U=0.D0 87 : DO 42 L=1,NPG 88 : U=U+FN(I,L)*FN(J,L)/RPG(L)/RPG(L)*PGSQ(L) 89 : 42 CONTINUE 90 : U=U*COEF(IC) 91 : AF2(KE,J,I)=AF2(KE,J,I)+U*AIMPL 92 : AS2(KE,J,I)=AS2(KE,J,I)+U*(AIMPL-1.D0) 93 : 41 CONTINUE 94 : ELSEIF(IAXI.EQ.2)THEN 95 : DO 43 I=1,NP 96 : DO 43 J=1,NP 97 : U=0.D0 98 : DO 44 L=1,NPG 99 : U=U+FN(I,L)*FN(J,L)/RPG(L)/RPG(L)*PGSQ(L) 100 : 44 CONTINUE 101 : U=U*COEF(IC) 102 : AF1(KE,J,I)=AF1(KE,J,I)+U*AIMPL 103 : AS1(KE,J,I)=AS1(KE,J,I)+U*(AIMPL-1.D0) 104 : 43 CONTINUE 105 : ENDIF 106 : ENDIF 107 : GO TO 108 108 : 109 : 5 CONTINUE 110 : DO 51 L=1,NPG 111 : C=0.D0 112 : DO 52 I=1,NP 113 : IU=LE(I,KE) 114 : C=C+COEF(IU)*FN(I,L) 115 : 52 CONTINUE 116 : CL(L)=C 117 : 51 CONTINUE 118 : C 119 : C write(6,*)' BCL 15 ' 120 : DO 15 I=1,NP 121 : DO 15 J=I,NP 122 : U=0.D0 123 : DO 12 L=1,NPG 124 : V=0.D0 125 : DO 13 N=1,IDIM 126 : V=V+HR(N,I,L)*HR(N,J,L) 127 : 13 CONTINUE 128 : U=U+V*PGSQ(L)*CL(L) 129 : 12 CONTINUE 130 : DO 14 K=1,NINC 131 : AF(J,I,K,K)=AF(J,I,K,K)+U 132 : IF(I.NE.J)AF(I,J,K,K)=AF(I,J,K,K)+U 133 : 14 CONTINUE 134 : IF(IAX.EQ.0)GO TO 15 135 : U=0.D0 136 : DO 410 L=1,NPG 137 : U=U+FN(I,L)*FN(J,L)/RPG(L)/RPG(L)*PGSQ(L)*CL(L) 138 : 410 CONTINUE 139 : KU=3-IAX 140 : AF(J,I,KU,KU)=AF(J,I,KU,KU)+U 141 : IF(I.NE.J)AF(I,J,KU,KU)=AF(I,J,KU,KU)+U 142 : 15 CONTINUE 143 : GO TO 108 144 : C 145 : C 146 : 70 CONTINUE 146 : IF(NINC.NE.1)CALL ARRET(0) 146 : IF(KJTT.EQ.4)CALL ARRET(0) 147 : IC=1+(1-IK1)*(NK-1) 148 : DO 71 I=1,NP 149 : DO 71 J=1,NP 150 : U=0.D0 151 : DO 72 L=1,NPG 152 : UL=0.D0 153 : DO 73 M=1,IDIM 154 : UN=0.D0 155 : DO 74 N=1,IDIM 156 : UN=UN+COES(N,M,IC)*HR(N,J,L) 157 : 74 CONTINUE 158 : UL=UL+UN*HR(M,I,L) 159 : 73 CONTINUE 160 : U=U+UL*PGSQ(L) 161 : 72 CONTINUE 162 : AF(J,I,1,1)=AF(J,I,1,1)+U 163 : 71 CONTINUE 164 : GO TO 108 165 : C 166 : C 167 : 80 CONTINUE 167 : IF(KJTT.EQ.4)CALL ARRET(0) 167 : IF(IARG.EQ.0)CALL ARRET(0) 168 : IC=1+(1-IK1)*(NK-1) 169 : ALT=COG(1,IC)-COG(2,IC) 170 : AT=COG(2,IC) 171 : C WRITE(6,*)' ALT ',ALT,' AT ',AT 172 : IF(IDIM.EQ.2) THEN 173 : CO=SQRT(S2(1,IC)*S2(1,IC)+S2(2,IC)*S2(2,IC)) 174 : ENDIF 175 : IF(IDIM.EQ.3) THEN 176 : CO=SQRT(S2(1,IC)*S2(1,IC)+S2(2,IC)*S2(2,IC)+S2(3,IC)*S2(3,IC)) 177 : ENDIF 178 : C WRITE(6,*)' CO ',CO 179 : C WRITE(6,*)' S2 ',S2(1,IC),S2(2,IC) 180 : IF(IDIM.EQ.2) THEN 181 : COE(1,1)=ALT*S2(1,IC)*S2(1,IC)/CO+AT*CO 182 : COE(2,2)=ALT*S2(2,IC)*S2(2,IC)/CO+AT*CO 183 : COE(1,2)=ALT*S2(1,IC)*S2(2,IC)/CO 184 : COE(2,1)=COE(1,2) 185 : ENDIF 186 : IF(IDIM.EQ.3) THEN 187 : COE(1,1)=ALT*S2(1,IC)*S2(1,IC)/CO+AT*CO 188 : COE(2,2)=ALT*S2(2,IC)*S2(2,IC)/CO+AT*CO 189 : COE(3,3)=ALT*S2(3,IC)*S2(3,IC)/CO+AT*CO 190 : COE(1,2)=ALT*S2(1,IC)*S2(2,IC)/CO 191 : COE(1,3)=ALT*S2(1,IC)*S2(3,IC)/CO 192 : COE(2,3)=ALT*S2(2,IC)*S2(3,IC)/CO 193 : COE(2,1)=COE(1,2) 194 : COE(3,1)=COE(1,3) 195 : COE(3,2)=COE(2,3) 196 : ENDIF 197 : C WRITE(6,*)' COE ',COE(1,1),COE(1,2),COE(1,3) 198 : C WRITE(6,*)' ',COE(2,1),COE(2,2),COE(2,3) 199 : C WRITE(6,*)' ',COE(3,1),COE(3,2),COE(3,3) 200 : C 201 : DO 81 I=1,NP 202 : DO 81 J=1,NP 203 : U=0.D0 204 : DO 82 L=1,NPG 205 : UL=0.D0 206 : DO 83 M=1,IDIM 207 : UN=0.D0 208 : DO 84 N=1,IDIM 209 : UN=UN+COE(N,M)*HR(N,J,L) 210 : 84 CONTINUE 211 : UL=UL+UN*HR(M,I,L) 212 : 83 CONTINUE 213 : U=U+UL*PGSQ(L) 214 : 82 CONTINUE 215 : AF(J,I,1,1)=AF(J,I,1,1)+U 216 : 81 CONTINUE 217 : 218 : 108 CONTINUE 219 : 220 : C write(6,*)' FIN XLAPL' 221 : 222 : RETURN 223 : 1002 FORMAT(10(1X,1PE11.4)) 224 : 1001 FORMAT(20(1X,I5)) 225 : END 226 : 227 : 228 : 229 : 230 : 231 : 232 :
© Cast3M 2003 - All rights reserved.
Disclaimer