Télécharger qsijs2.eso

Retour à la liste

Numérotation des lignes :

  1. C QSIJS2 SOURCE BP208322 16/11/18 21:20:37 9177
  2. SUBROUTINE QSIJS2(XEL,MELE,NBMA1,IDI,XPI,SHP,QSI,IROUT)
  3. C======================================================================C
  4. C C
  5. C INPUT: C
  6. C XEL(3,*) : LES COORDONNEES DE L'ELEMENT DANS LE REPERE LOCAL C
  7. C MELE : NUMERO DE L ELEMENT C
  8. C NBMA1 NOMBRE DE NOEUDS DE L ELEMENT C
  9. C XPI(3): COORDONNEES DU POINT DANS LE SYSTEME GLOBAL C
  10. C IDI 2 POUR 2D ET ELEMENTS SURFACIQUES 3 POUR ELEM. VOLUMIQUES C
  11. C C
  12. C INPUT/OUTPUT: C
  13. C QSI(3): COORDONNEES DU POINT DANS L ELEMENT DE REFERENCE C
  14. C C
  15. C OUTPUT: C
  16. C C
  17. C SHP FONCTIONS DE FORMES ET DERIVEES C
  18. C IROUT 0 SI SUCCES DE L OPERATION 1 SI ECHEC C
  19. C Sortie convergence a 1.E-5 sur élement de référence C
  20. C J.S. FLEURET / TRANSFORME JM BAZE : POUR TOUS ELEMENTS(EN PRINCIPE) C
  21. C======================================================================C
  22. C Gounand 13/07/2007 : la différence par rapport à qsijs est qu'on
  23. C fournit une estimation initiale du tableau QSI
  24. IMPLICIT INTEGER(I-N)
  25. IMPLICIT REAL*8 (A-H,O-Z)
  26. -INC CCOPTIO
  27. -INC CCGEOME
  28. DIMENSION XEL(3,*),QSI(3),XPI(3),DQS(3)
  29. DIMENSION XSU(20),YSU(20),ZSU(20)
  30. DIMENSION SHP(6,*),D(3,3),DINV(3,3),V(3),XU(3),F(3)
  31. EXTERNAL SHAPE
  32. DATA UN,XZER,EPSI/1.D0,0.D0,1.D-7/
  33.  
  34. C- Quelques initialisations
  35. IROUT=0
  36. Iter=0
  37. IAcc=0
  38. IA1=4
  39.  
  40. * QSI(1)=XZER
  41. * QSI(2)=XZER
  42. * QSI(3)=XZER
  43.  
  44. IF (IDI.EQ.1) GOTO 800
  45.  
  46. DXDQSI=XZER
  47. DXDETA=XZER
  48. DYDQSI=XZER
  49. DYDETA=XZER
  50.  
  51. C-------------------- DEBUT DES ITERATIONS ------------------------
  52. 600 CONTINUE
  53. ITER = ITER + 1
  54. IACC = IACC + 1
  55. IF (ITER.GT.100) THEN
  56. IROUT = ITER
  57. RETURN
  58. ENDIF
  59. * WRITE(IOIMP,*) ' qsijs2 : QSI1=',QSI(1),' QSI2=',QSI(2)
  60. CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,IRET)
  61. if (iret.eq.0) then
  62. MOTERR(1:4)=NOMS(MELE)
  63. CALL ERREUR(68)
  64. RETURN
  65. endif
  66.  
  67.  
  68. C PROJECTION
  69. DO 501 IK = 1,IDI
  70. XU(IK)= XZER
  71. DO 501 JK = 1,NBMA1
  72. XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK))
  73. 501 CONTINUE
  74.  
  75. C DD = XZER
  76. DO 502 IK = 1,IDI
  77. F(IK) = XU(IK) - XPI(IK)
  78. 502 CONTINUE
  79. C
  80. C CAS 2 ELEMENTS SURFACIQUES
  81. C
  82. IF (IDI.EQ.2) THEN
  83. C IF ((MELE.LT.14).OR.(MELE.GT.26)) THEN
  84. DO 100 I=1,NBMA1
  85. DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I)
  86. DXDETA=DXDETA+SHP(3,I)*XEL(1,I)
  87. DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I)
  88. DYDETA=DYDETA+SHP(3,I)*XEL(2,I)
  89. 100 CONTINUE
  90. DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI
  91. XXXX = DJAC
  92. IF(DJAC.NE.XZER) XXXX=UN/DJAC
  93. DQSIDX= DYDETA*XXXX
  94. DQSIDY=-DXDETA*XXXX
  95. DETADX=-DYDQSI*XXXX
  96. DETADY= DXDQSI*XXXX
  97. DINV(1,1) = DQSIDX
  98. DINV(2,1) = DQSIDY
  99. DINV(1,2) = DETADX
  100. DINV(2,2) = DETADY
  101.  
  102. ELSE
  103. C
  104. C ELEMENTS VOLUMIQUES
  105. C
  106.  
  107. DO 310 I=1,3
  108. DO 310 J=1,3
  109. D(I,J)=XZER
  110. 310 CONTINUE
  111. DO 300 I=1,NBMA1
  112. C
  113. D(1,1)=D(1,1)+SHP(2,I)*XEL(1,I)
  114. D(2,1)=D(2,1)+SHP(3,I)*XEL(1,I)
  115. D(3,1)=D(3,1)+SHP(4,I)*XEL(1,I)
  116. C
  117. D(1,2)=D(1,2)+SHP(2,I)*XEL(2,I)
  118. D(2,2)=D(2,2)+SHP(3,I)*XEL(2,I)
  119. D(3,2)=D(3,2)+SHP(4,I)*XEL(2,I)
  120. C
  121. D(1,3)=D(1,3)+SHP(2,I)*XEL(3,I)
  122. D(2,3)=D(2,3)+SHP(3,I)*XEL(3,I)
  123. D(3,3)=D(3,3)+SHP(4,I)*XEL(3,I)
  124. 300 CONTINUE
  125. C INVERSION DE LA MATRICE D(3,3)
  126. DINV(1,1)= D(2,2)*D(3,3)-D(2,3)*D(3,2)
  127. DINV(2,2)= D(1,1)*D(3,3)-D(1,3)*D(3,1)
  128. DINV(3,3)= D(1,1)*D(2,2)-D(1,2)*D(2,1)
  129. DINV(1,2)=-D(1,2)*D(3,3)+D(3,2)*D(1,3)
  130. DINV(2,1)=-D(2,1)*D(3,3)+D(2,3)*D(3,1)
  131. DINV(1,3)= D(1,2)*D(2,3)-D(2,2)*D(1,3)
  132. DINV(3,1)= D(2,1)*D(3,2)-D(2,2)*D(3,1)
  133. DINV(2,3)=-D(1,1)*D(2,3)+D(2,1)*D(1,3)
  134. DINV(3,2)=-D(1,1)*D(3,2)+D(1,2)*D(3,1)
  135. DJAC=D(1,1)*DINV(1,1)+D(2,1)*DINV(1,2)+D(3,1)*DINV(1,3)
  136. XXXX = DJAC
  137. IF(DJAC.NE.XZER) XXXX=UN/DJAC
  138. DO 350 IA=1,3
  139. DO 350 IB=1,3
  140. DINV(IA,IB)=DINV(IA,IB)*XXXX
  141. 350 CONTINUE
  142. ENDIF
  143. C
  144. DDQ =XZER
  145. DO 505 I=1,IDI
  146. DQS(I)=XZER
  147. DO 506 J=1,IDI
  148. DQS(I) =DQS(I)+ DINV(J,I)*F(J)
  149. 506 CONTINUE
  150. DDQ = DDQ + DQS(I)*DQS(I)
  151. 505 QSI(I) = QSI(I)-DQS(I)
  152. DDQ = SQRT(DDQ)
  153. IF(DDQ.LT.EPSI) GOTO 699
  154. XSU(IACC) = QSI(1)
  155. YSU(IACC) = QSI(2)
  156. IF(IDI.EQ.3) ZSU(IACC) = QSI(3)
  157. C
  158. IF(IACC.GT.IA1) THEN
  159. C acceleration
  160. dqs1 = 2.D0*XSU(IACC-1)-QSI(1)-XSU(IACC-2)
  161. if (abs(dqs1).GT.1.D-10) then
  162. BX = (QSI(1) - XSU(IACC-1))*2.D0/dqs1
  163. QSI(1) =(1.D0+BX)*XSU(IACC-1) - BX*XSU(IACC-2)
  164. end if
  165. dqs2 = 2.D0*YSU(IACC-1)-QSI(2)-YSU(IACC-2)
  166. if (abs(dqs2).GT.1.D-10) then
  167. BY = (QSI(2) - YSU(IACC-1))*2.D0/dqs2
  168. QSI(2) =(1.D0+BY)*YSU(IACC-1) - BY*YSU(IACC-2)
  169. endif
  170. XSU(IACC) = QSI(1)
  171. YSU(IACC) = QSI(2)
  172. IF(IDI.EQ.3) THEN
  173. dqs3 = 2.D0*ZSU(IACC-1)-QSI(3)-ZSU(IACC-2)
  174. if (abs(dqs3).GT.1.D-10) then
  175. BZ = (QSI(3) - ZSU(IACC-1))*2.D0/dqs3
  176. QSI(3) =(1.D0+BZ)*ZSU(IACC-1) - BZ*ZSU(IACC-2)
  177. endif
  178. ZSU(IACC) = QSI(3)
  179. ENDIF
  180. IA1 = IA1+2
  181. ENDIF
  182. C au cas ou il faudrait plus de 20 acceleration
  183. if ( iacc.EQ.20) then
  184. do 400 i = 1,2
  185. xsu(i) = xsu(iacc-2+i)
  186. ysu(i) = ysu(iacc-2+i)
  187. zsu(i) = zsu(iacc-2+i)
  188. 400 continue
  189. ia1 = 1
  190. iacc =0
  191. endif
  192.  
  193. C
  194. C-------- ITERATIONS ----------------------
  195.  
  196. GOTO 600
  197. * ENDIF
  198.  
  199. C========== Cas IDI = 1 ==========
  200. C Methode de Newton-Raphson pour resoudre X(QSI)-XPI(1)=0
  201. 800 Iter=Iter+1
  202. IAcc=IAcc+1
  203. IF (Iter.GT.100) THEN
  204. IROUT=Iter
  205. RETURN
  206. ENDIF
  207. CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,IRET)
  208. if (iret.eq.0) then
  209. MOTERR(1:4)=NOMS(MELE)
  210. CALL ERREUR(68)
  211. RETURN
  212. endif
  213. XU1=XZER
  214. dXdQSI=XZER
  215. DO i=1,NBMA1
  216. XU1=XU1+SHP(1,i)*XEL(1,i)
  217. dXdQSI=dXdQSI+SHP(2,i)*XEL(1,i)
  218. ENDDO
  219. IF (dXdQSI.NE.XZER) THEN
  220. F1=XU1-XPI(1)
  221. dQS1=F1/dXdQSI
  222. QSI(1)=QSI(1)-dQS1
  223. ELSE
  224. dQS1=XZER
  225. ENDIF
  226. C Test de convergence
  227. IF (ABS(dQS1).LT.EPSI) RETURN
  228. C Acceleration
  229. IF (IAcc.GT.IA1) THEN
  230. dQS1=2.*XSU(IAcc-1)-QSI(1)-XSU(IAcc-2)
  231. IF (ABS(dQS1).GT.1.D-10) THEN
  232. BX=2.*(QSI(1)-XSU(IAcc-1))/dQS1
  233. QSI(1)=(UN+BX)*XSU(IAcc-1)-BX*XSU(IAcc-2)
  234. ENDIF
  235. IA1=IA1+2
  236. ENDIF
  237. XSU(IAcc)=QSI(1)
  238. C Dans le cas ou il faudrait plus de 20 accelerations
  239. IF (IAcc.EQ.20) THEN
  240. XSU(1)=XSU(IAcc-1)
  241. XSU(2)=XSU(IAcc)
  242. IA1=1
  243. IAcc=0
  244. ENDIF
  245. GOTO 800
  246.  
  247. 699 CONTINUE
  248. RETURN
  249. END
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  

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