Télécharger qsijs.eso

Retour à la liste

Numérotation des lignes :

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

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