Télécharger qsijs2.eso

Retour à la liste

Numérotation des lignes :

  1. C QSIJS2 SOURCE CB215821 19/03/18 21:16:10 10160
  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.  
  27. -INC PPARAM
  28. -INC CCOPTIO
  29. -INC CCGEOME
  30. DIMENSION XEL(3,*),QSI(3),XPI(3),DQS(3)
  31. DIMENSION XSU(20),YSU(20),ZSU(20)
  32. DIMENSION SHP(6,*),D(3,3),DINV(3,3),V(3),XU(3),F(3)
  33. EXTERNAL SHAPE
  34. DATA UN,XZER,EPSI/1.D0,0.D0,1.D-7/
  35.  
  36. C- Quelques initialisations
  37. IROUT=0
  38. Iter=0
  39. IAcc=0
  40. IA1=4
  41.  
  42. * QSI(1)=XZER
  43. * QSI(2)=XZER
  44. * QSI(3)=XZER
  45.  
  46. IF (IDI.EQ.1) GOTO 800
  47.  
  48. DXDQSI=XZER
  49. DXDETA=XZER
  50. DYDQSI=XZER
  51. DYDETA=XZER
  52.  
  53. C-------------------- DEBUT DES ITERATIONS ------------------------
  54. 600 CONTINUE
  55. ITER = ITER + 1
  56. IACC = IACC + 1
  57. IF (ITER.GT.100) THEN
  58. IROUT = ITER
  59. RETURN
  60. ENDIF
  61. * WRITE(IOIMP,*) ' qsijs2 : QSI1=',QSI(1),' QSI2=',QSI(2)
  62. CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,IRET)
  63. if (iret.eq.0) then
  64. MOTERR(1:4)=NOMS(MELE)
  65. CALL ERREUR(68)
  66. RETURN
  67. endif
  68.  
  69.  
  70. C PROJECTION
  71. DO 501 IK = 1,IDI
  72. XU(IK)= XZER
  73. DO 501 JK = 1,NBMA1
  74. XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK))
  75. 501 CONTINUE
  76.  
  77. C DD = XZER
  78. DO 502 IK = 1,IDI
  79. F(IK) = XU(IK) - XPI(IK)
  80. 502 CONTINUE
  81. C
  82. C CAS 2 ELEMENTS SURFACIQUES
  83. C
  84. IF (IDI.EQ.2) THEN
  85. C IF ((MELE.LT.14).OR.(MELE.GT.26)) THEN
  86. DO 100 I=1,NBMA1
  87. DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I)
  88. DXDETA=DXDETA+SHP(3,I)*XEL(1,I)
  89. DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I)
  90. DYDETA=DYDETA+SHP(3,I)*XEL(2,I)
  91. 100 CONTINUE
  92. DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI
  93. XXXX = DJAC
  94. IF(DJAC.NE.XZER) XXXX=UN/DJAC
  95. DQSIDX= DYDETA*XXXX
  96. DQSIDY=-DXDETA*XXXX
  97. DETADX=-DYDQSI*XXXX
  98. DETADY= DXDQSI*XXXX
  99. DINV(1,1) = DQSIDX
  100. DINV(2,1) = DQSIDY
  101. DINV(1,2) = DETADX
  102. DINV(2,2) = DETADY
  103.  
  104. ELSE
  105. C
  106. C ELEMENTS VOLUMIQUES
  107. C
  108.  
  109. DO 310 I=1,3
  110. DO 310 J=1,3
  111. D(I,J)=XZER
  112. 310 CONTINUE
  113. DO 300 I=1,NBMA1
  114. C
  115. D(1,1)=D(1,1)+SHP(2,I)*XEL(1,I)
  116. D(2,1)=D(2,1)+SHP(3,I)*XEL(1,I)
  117. D(3,1)=D(3,1)+SHP(4,I)*XEL(1,I)
  118. C
  119. D(1,2)=D(1,2)+SHP(2,I)*XEL(2,I)
  120. D(2,2)=D(2,2)+SHP(3,I)*XEL(2,I)
  121. D(3,2)=D(3,2)+SHP(4,I)*XEL(2,I)
  122. C
  123. D(1,3)=D(1,3)+SHP(2,I)*XEL(3,I)
  124. D(2,3)=D(2,3)+SHP(3,I)*XEL(3,I)
  125. D(3,3)=D(3,3)+SHP(4,I)*XEL(3,I)
  126. 300 CONTINUE
  127. C INVERSION DE LA MATRICE D(3,3)
  128. DINV(1,1)= D(2,2)*D(3,3)-D(2,3)*D(3,2)
  129. DINV(2,2)= D(1,1)*D(3,3)-D(1,3)*D(3,1)
  130. DINV(3,3)= D(1,1)*D(2,2)-D(1,2)*D(2,1)
  131. DINV(1,2)=-D(1,2)*D(3,3)+D(3,2)*D(1,3)
  132. DINV(2,1)=-D(2,1)*D(3,3)+D(2,3)*D(3,1)
  133. DINV(1,3)= D(1,2)*D(2,3)-D(2,2)*D(1,3)
  134. DINV(3,1)= D(2,1)*D(3,2)-D(2,2)*D(3,1)
  135. DINV(2,3)=-D(1,1)*D(2,3)+D(2,1)*D(1,3)
  136. DINV(3,2)=-D(1,1)*D(3,2)+D(1,2)*D(3,1)
  137. DJAC=D(1,1)*DINV(1,1)+D(2,1)*DINV(1,2)+D(3,1)*DINV(1,3)
  138. XXXX = DJAC
  139. IF(DJAC.NE.XZER) XXXX=UN/DJAC
  140. DO 350 IA=1,3
  141. DO 350 IB=1,3
  142. DINV(IA,IB)=DINV(IA,IB)*XXXX
  143. 350 CONTINUE
  144. ENDIF
  145. C
  146. DDQ =XZER
  147. DO 505 I=1,IDI
  148. DQS(I)=XZER
  149. DO 506 J=1,IDI
  150. DQS(I) =DQS(I)+ DINV(J,I)*F(J)
  151. 506 CONTINUE
  152. DDQ = DDQ + DQS(I)*DQS(I)
  153. 505 QSI(I) = QSI(I)-DQS(I)
  154. DDQ = SQRT(DDQ)
  155. IF(DDQ.LT.EPSI) GOTO 699
  156. XSU(IACC) = QSI(1)
  157. YSU(IACC) = QSI(2)
  158. IF(IDI.EQ.3) ZSU(IACC) = QSI(3)
  159. C
  160. IF(IACC.GT.IA1 .AND. IACC.GT.2) THEN
  161. C acceleration
  162. dqs1 = 2.D0*XSU(IACC-1)-QSI(1)-XSU(IACC-2)
  163. if (abs(dqs1).GT.1.D-10) then
  164. BX = (QSI(1) - XSU(IACC-1))*2.D0/dqs1
  165. QSI(1) =(1.D0+BX)*XSU(IACC-1) - BX*XSU(IACC-2)
  166. end if
  167. dqs2 = 2.D0*YSU(IACC-1)-QSI(2)-YSU(IACC-2)
  168. if (abs(dqs2).GT.1.D-10) then
  169. BY = (QSI(2) - YSU(IACC-1))*2.D0/dqs2
  170. QSI(2) =(1.D0+BY)*YSU(IACC-1) - BY*YSU(IACC-2)
  171. endif
  172. XSU(IACC) = QSI(1)
  173. YSU(IACC) = QSI(2)
  174. IF(IDI.EQ.3) THEN
  175. dqs3 = 2.D0*ZSU(IACC-1)-QSI(3)-ZSU(IACC-2)
  176. if (abs(dqs3).GT.1.D-10) then
  177. BZ = (QSI(3) - ZSU(IACC-1))*2.D0/dqs3
  178. QSI(3) =(1.D0+BZ)*ZSU(IACC-1) - BZ*ZSU(IACC-2)
  179. endif
  180. ZSU(IACC) = QSI(3)
  181. ENDIF
  182. IA1 = IA1+2
  183. ENDIF
  184. C au cas ou il faudrait plus de 20 acceleration
  185. if ( iacc.EQ.20) then
  186. do 400 i = 1,2
  187. xsu(i) = xsu(iacc-2+i)
  188. ysu(i) = ysu(iacc-2+i)
  189. zsu(i) = zsu(iacc-2+i)
  190. 400 continue
  191. ia1 = 1
  192. iacc =0
  193. endif
  194.  
  195. C
  196. C-------- ITERATIONS ----------------------
  197.  
  198. GOTO 600
  199. * ENDIF
  200.  
  201. C========== Cas IDI = 1 ==========
  202. C Methode de Newton-Raphson pour resoudre X(QSI)-XPI(1)=0
  203. 800 Iter=Iter+1
  204. IAcc=IAcc+1
  205. IF (Iter.GT.100) THEN
  206. IROUT=Iter
  207. RETURN
  208. ENDIF
  209. CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,IRET)
  210. if (iret.eq.0) then
  211. MOTERR(1:4)=NOMS(MELE)
  212. CALL ERREUR(68)
  213. RETURN
  214. endif
  215. XU1=XZER
  216. dXdQSI=XZER
  217. DO i=1,NBMA1
  218. XU1=XU1+SHP(1,i)*XEL(1,i)
  219. dXdQSI=dXdQSI+SHP(2,i)*XEL(1,i)
  220. ENDDO
  221. IF (dXdQSI.NE.XZER) THEN
  222. F1=XU1-XPI(1)
  223. dQS1=F1/dXdQSI
  224. QSI(1)=QSI(1)-dQS1
  225. ELSE
  226. dQS1=XZER
  227. ENDIF
  228. C Test de convergence
  229. IF (ABS(dQS1).LT.EPSI) RETURN
  230. C Acceleration
  231. IF (IAcc.GT.IA1) THEN
  232. dQS1=2.*XSU(IAcc-1)-QSI(1)-XSU(IAcc-2)
  233. IF (ABS(dQS1).GT.1.D-10) THEN
  234. BX=2.*(QSI(1)-XSU(IAcc-1))/dQS1
  235. QSI(1)=(UN+BX)*XSU(IAcc-1)-BX*XSU(IAcc-2)
  236. ENDIF
  237. IA1=IA1+2
  238. ENDIF
  239. XSU(IAcc)=QSI(1)
  240. C Dans le cas ou il faudrait plus de 20 accelerations
  241. IF (IAcc.EQ.20) THEN
  242. XSU(1)=XSU(IAcc-1)
  243. XSU(2)=XSU(IAcc)
  244. IA1=1
  245. IAcc=0
  246. ENDIF
  247. GOTO 800
  248.  
  249. 699 CONTINUE
  250. END
  251.  
  252.  

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