Télécharger kdom7.eso

Retour à la liste

Numérotation des lignes :

  1. C KDOM7 SOURCE KK2000 14/04/10 21:15:14 8032
  2. SUBROUTINE KDOM7(IP1,IP2,IP3,IP4,XCEN)
  3. C
  4. C************************************************************************
  5. C
  6. C PROJET : CASTEM 2000
  7. C
  8. C NOM : KDOM7
  9. C
  10. C DESCRIPTION : Subroutine called by LEKMA0
  11. C Given a QUA4 (IP1,IP2,IP3,IP4), this subroutine
  12. C compute its center by averaging the centers of
  13. C the two possible configuration it can be
  14. C divided in
  15. C
  16. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  17. C
  18. C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF
  19. C
  20. C************************************************************************
  21. C
  22. C
  23.  
  24. IMPLICIT INTEGER(I-N)
  25. IMPLICIT REAL*8(A-H,O-Z)
  26.  
  27. INTEGER IP1,IP2,IP3, IP4, I1,I2, IPCEN, IPL, IPR
  28. REAL*8 VOL, XCEN(3),VOLT,DX0,DY0,DZ0,DX1,DY1,DZ1,CEL
  29. & ,RNOX0,RNOY0,RNOZ0,RNOX1,RNOY1,RNOZ1,XCEN1(3),VOLT1,PSCA
  30. & ,RNOT
  31. LOGICAL LOGREG
  32. -INC CCOPTIO
  33. -INC SMCOORD
  34. C
  35. DZ0=0.0D0
  36. DZ1=0.0D0
  37. RNOT=0.0D0
  38. C
  39. C**** We chose as positive direction the one which
  40. C has the maximum modulus IP1IP(1+1) x IP1IP(1+3)
  41. C
  42. DO I1=1,4
  43. IF(I1 .EQ. 1)THEN
  44. C We consider the triangles IP1, IP2, IP4
  45. IPCEN=IP1
  46. IPR=IP2
  47. IPL=IP4
  48. ELSEIF(I1 .EQ. 2)THEN
  49. IPCEN=IP2
  50. IPR=IP3
  51. IPL=IP1
  52. ELSEIF(I1 .EQ. 3)THEN
  53. IPCEN=IP3
  54. IPR=IP4
  55. IPL=IP2
  56. ELSEIF(I1 .EQ. 4)THEN
  57. IPCEN=IP4
  58. IPR=IP1
  59. IPL=IP3
  60. ENDIF
  61. DX0=XCOOR((IPR-1)*(IDIM+1)+1)-
  62. & XCOOR((IPCEN-1)*(IDIM+1)+1)
  63. DY0=XCOOR((IPR-1)*(IDIM+1)+2)-
  64. & XCOOR((IPCEN-1)*(IDIM+1)+2)
  65. IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-
  66. & XCOOR((IPCEN-1)*(IDIM+1)+3)
  67. DX1=XCOOR((IPL-1)*(IDIM+1)+1)-
  68. & XCOOR((IPCEN-1)*(IDIM+1)+1)
  69. DY1=XCOOR((IPL-1)*(IDIM+1)+2)-
  70. & XCOOR((IPCEN-1)*(IDIM+1)+2)
  71. IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-
  72. & XCOOR((IPCEN-1)*(IDIM+1)+3)
  73. RNOX1=(DY0*DZ1-DY1*DZ0)
  74. RNOY1=(DZ0*DX1-DZ1*DX0)
  75. RNOZ1=(DX0*DY1-DX1*DY0)
  76. VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0
  77. VOL=VOL/0.5D0
  78. IF(VOL .GT. RNOT)THEN
  79. RNOT=VOL
  80. RNOX0=RNOX1
  81. RNOY0=RNOY1
  82. RNOZ0=RNOZ1
  83. ENDIF
  84. ENDDO
  85. C
  86. C**** First division in triangles
  87. C We consider the triangles IP3, IP1, IP2 and
  88. C IP4, IP1, IP3
  89. C
  90. XCEN1(1)=0.0D0
  91. XCEN1(2)=0.0D0
  92. XCEN1(3)=0.0D0
  93. XCEN(1)=0.0D0
  94. XCEN(2)=0.0D0
  95. XCEN(3)=0.0D0
  96. VOLT=0.0D0
  97. VOLT1=0.0D0
  98. LOGREG=.TRUE.
  99. C
  100. DO I1=1,2
  101. IF(I1 .EQ. 1)THEN
  102. C We consider the triangles IP3, IP1, IP2
  103. IPCEN=IP1
  104. IPR=IP2
  105. IPL=IP3
  106. ELSEIF(I1 .EQ. 2)THEN
  107. IPCEN=IP1
  108. IPR=IP3
  109. IPL=IP4
  110. ENDIF
  111. DX0=XCOOR((IPR-1)*(IDIM+1)+1)-
  112. & XCOOR((IPCEN-1)*(IDIM+1)+1)
  113. DY0=XCOOR((IPR-1)*(IDIM+1)+2)-
  114. & XCOOR((IPCEN-1)*(IDIM+1)+2)
  115. IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-
  116. & XCOOR((IPCEN-1)*(IDIM+1)+3)
  117. DX1=XCOOR((IPL-1)*(IDIM+1)+1)-
  118. & XCOOR((IPCEN-1)*(IDIM+1)+1)
  119. DY1=XCOOR((IPL-1)*(IDIM+1)+2)-
  120. & XCOOR((IPCEN-1)*(IDIM+1)+2)
  121. IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-
  122. & XCOOR((IPCEN-1)*(IDIM+1)+3)
  123. RNOX1=(DY0*DZ1-DY1*DZ0)
  124. RNOY1=(DZ0*DX1-DZ1*DX0)
  125. RNOZ1=(DX0*DY1-DX1*DY0)
  126. VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0
  127. VOL=VOL*0.5D0
  128. PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1)
  129. IF((VOL .LE. (RNOT*1.0D-6)) .OR. (PSCA.LE.0))THEN
  130. LOGREG=.FALSE.
  131. ENDIF
  132. VOLT1=VOLT1+VOL
  133. DO I2=1,IDIM,1
  134. CEL=(XCOOR((IPCEN-1)*(IDIM+1)+I2)+
  135. & XCOOR((IPL-1)*(IDIM+1)+I2)+XCOOR((IPR-1)*(IDIM+1)+I2))
  136. & /3.0D0
  137. XCEN1(I2)=XCEN1(I2)+(VOL*CEL)
  138. ENDDO
  139. ENDDO
  140. C
  141. C**** If one of the area is negative or almost zero the division is not
  142. C possible (LOGREG=.FALSE.)
  143. C
  144. IF(LOGREG)THEN
  145. VOLT=VOLT1
  146. DO I2=1,IDIM,1
  147. XCEN(I2)=XCEN1(I2)
  148. ENDDO
  149. ENDIF
  150. C
  151. C**** Second division in triangles
  152. C We consider the triangles IP3, IP2, IP4 and
  153. C IP4, IP2, IP1
  154. C
  155. C
  156. VOLT1=0.0D0
  157. LOGREG=.TRUE.
  158. XCEN1(1)=0.0D0
  159. XCEN1(2)=0.0D0
  160. XCEN1(3)=0.0D0
  161. C
  162. DO I1=1,2
  163. IF(I1 .EQ. 1)THEN
  164. C We consider the triangles
  165. IPCEN=IP2
  166. IPR=IP3
  167. IPL=IP4
  168. ELSEIF(I1 .EQ. 2)THEN
  169. IPCEN=IP2
  170. IPR=IP4
  171. IPL=IP1
  172. ENDIF
  173. DX0=XCOOR((IPR-1)*(IDIM+1)+1)-
  174. & XCOOR((IPCEN-1)*(IDIM+1)+1)
  175. DY0=XCOOR((IPR-1)*(IDIM+1)+2)-
  176. & XCOOR((IPCEN-1)*(IDIM+1)+2)
  177. IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-
  178. & XCOOR((IPCEN-1)*(IDIM+1)+3)
  179. DX1=XCOOR((IPL-1)*(IDIM+1)+1)-
  180. & XCOOR((IPCEN-1)*(IDIM+1)+1)
  181. DY1=XCOOR((IPL-1)*(IDIM+1)+2)-
  182. & XCOOR((IPCEN-1)*(IDIM+1)+2)
  183. IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-
  184. & XCOOR((IPCEN-1)*(IDIM+1)+3)
  185. RNOX1=(DY0*DZ1-DY1*DZ0)
  186. RNOY1=(DZ0*DX1-DZ1*DX0)
  187. RNOZ1=(DX0*DY1-DX1*DY0)
  188. VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0
  189. VOL=VOL*0.5D0
  190. PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1)
  191. IF((VOL .LE. (RNOT*1.0D-6)) .OR. (PSCA.LE.0))THEN
  192. LOGREG=.FALSE.
  193. ENDIF
  194. VOLT1=VOLT1+VOL
  195. DO I2=1,IDIM,1
  196. CEL=(XCOOR((IPCEN-1)*(IDIM+1)+I2)+
  197. & XCOOR((IPL-1)*(IDIM+1)+I2)+XCOOR((IPR-1)*(IDIM+1)+I2))
  198. & /3.0D0
  199. XCEN1(I2)=XCEN1(I2)+(VOL*CEL)
  200. ENDDO
  201. ENDDO
  202. C
  203. IF(LOGREG)THEN
  204. VOLT=VOLT+VOLT1
  205. DO I2=1,IDIM,1
  206. XCEN(I2)=XCEN(I2)+XCEN1(I2)
  207. ENDDO
  208. ENDIF
  209. C
  210. DO I2=1,IDIM,1
  211. XCEN(I2)=XCEN(I2)/VOLT
  212. ENDDO
  213. C
  214. C**** Finally, we check whether the center of gravity is "inside"
  215. C the domain, i.e. whether the normals CENIPi x CENIP(i+1)
  216. C have the same direction of (RNOX0,RNOY0,RNOZ0)
  217. C
  218. LOGREG=.FALSE.
  219. DO I1=1,4
  220. IF(I1 .EQ. 1)THEN
  221. C We consider the triangles CEN, IP1, IP2
  222. IPR=IP1
  223. IPL=IP2
  224. ELSEIF(I1 .EQ. 2)THEN
  225. IPR=IP2
  226. IPL=IP3
  227. ELSEIF(I1 .EQ. 3)THEN
  228. IPR=IP3
  229. IPL=IP4
  230. ELSEIF(I1 .EQ. 4)THEN
  231. IPR=IP4
  232. IPL=IP1
  233. ENDIF
  234. DX0=XCOOR((IPR-1)*(IDIM+1)+1)-
  235. & XCEN(1)
  236. DY0=XCOOR((IPR-1)*(IDIM+1)+2)-
  237. & XCEN(2)
  238. IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-
  239. & XCEN(3)
  240. DX1=XCOOR((IPL-1)*(IDIM+1)+1)-
  241. & XCEN(1)
  242. DY1=XCOOR((IPL-1)*(IDIM+1)+2)-
  243. & XCEN(2)
  244. IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-
  245. & XCEN(3)
  246. RNOX1=(DY0*DZ1-DY1*DZ0)
  247. RNOY1=(DZ0*DX1-DZ1*DX0)
  248. RNOZ1=(DX0*DY1-DX1*DY0)
  249. PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1)
  250. IF(PSCA .LE. 0.0D0)THEN
  251. LOGREG=.TRUE.
  252. ENDIF
  253. ENDDO
  254. C
  255. IF(LOGREG)THEN
  256. C
  257. WRITE(IOIMP,*) 'Subroutine kdom7.eso'
  258. WRITE(IOIMP,*) 'QUA4'
  259. C
  260. C Erreur -107 0
  261. C Liste des noeuds de l'élément :
  262. C
  263. CALL ERREUR(-107)
  264. WRITE(IOIMP,*) IP1, IP2, IP3, IP4
  265. C
  266. C Erreur 845 2
  267. C Maillage donné incorrect ?!!!
  268. C
  269. CALL ERREUR(845)
  270. C
  271. ENDIF
  272. C
  273. RETURN
  274. END
  275.  
  276.  
  277.  
  278.  

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