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

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