Télécharger qsijs.eso

Retour à la liste

Numérotation des lignes :

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

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