Télécharger melim.eso

Retour à la liste

Numérotation des lignes :

  1. C MELIM SOURCE PV 16/11/17 22:00:43 9180
  2. SUBROUTINE MELIM(MATRIK,KCLIM,
  3. $ INCX,KS2B,
  4. $ KMORS,KISA,
  5. $ IMPR,IRET)
  6. IMPLICIT REAL*8 (A-H,O-Z)
  7. IMPLICIT INTEGER (I-N)
  8. C***********************************************************************
  9. C NOM : MELIM
  10. C DESCRIPTION :
  11. C Ce SP s'occupe de prendre en compte les conditions
  12. C aux limites de type Dirichlet pour un système linéaire
  13. C du type Ax=b (avec A utilisant le stockage Morse.
  14. C On procède en trois étapes :
  15. C 1 - on surcharge le second membre (b)
  16. C et l'inconnue (x) avec le chpoint des conditions
  17. C aux limites
  18. C 2 - on s'occupe de la matrice Morse (A)
  19. C Aux degrés de libertés où il y a des
  20. C conditions aux limites : 1 sur la diagonale
  21. C 0 sur la ligne
  22. C 0 sur la colonne en reportant les
  23. C produits dans le second membre
  24. C (sinon perte de symétrie
  25. C éventuelle de A)
  26. C 3 - On nettoie la matrice morse A des 0 qu'on vient de lui
  27. C rajouter (sous-programme clmors).
  28. C
  29. C 2 et 3 peuvent sans doute etre regroupées au prix d'une
  30. C perte de lisibilité
  31. C
  32. C
  33. C LANGAGE : ESOPE
  34. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  35. C mél : gounand@semt2.smts.cea.fr
  36. C***********************************************************************
  37. C APPELES : KRIPAD, CLMORS
  38. C***********************************************************************
  39. C ENTREES : MATRIK, KCLIM, IMPR
  40. C ENTREES/SORTIES : KS2B, INCX
  41. C SORTIES : KMORS,KISA,IRET
  42. C CODE RETOUR (IRET) : 0 si ok
  43. C <0 si problème
  44. C MATRIK : pointeur sur segment MATRIK de l'include SMMATRIK
  45. C on pioche dedans les informations nécessaires
  46. C (différents pointeurs, nb. de ddl...)
  47. C KCLIM : pointeur sur segment MCHPOI de l'include SMCHPOI
  48. C chpoint contenant les conditions aux limites de
  49. C Dirichlet (valeur imposée).
  50. C IMPR : niveau d'impression
  51. C KS2B : pointeur sur segment IZA de l'include SMMATRIK
  52. C vecteur second membre
  53. C (ie b dans la résolution de Ax=b).
  54. C INCX : pointeur sur segment IZA de l'include SMMATRIK
  55. C vecteur des inconnues
  56. C (ie x dans la résolution de Ax=b).
  57. C KMORS : pointeurs sur segment PMORS et IZA de l'include
  58. C KISA SMMATRIK : profil et valeurs de la matrice morse
  59. C (ie A dans la résolution de Ax=b).
  60. C***********************************************************************
  61. C VERSION : v1, 20/12/99
  62. C HISTORIQUE : v1, 01/04/98, création
  63. C HISTORIQUE : 20/12/99,
  64. C Nouvelle signification de NUAN (renumérotation des ddl et non pas
  65. C des points)
  66. C HISTORIQUE : 05/01/00 : non-duplication de la matrice assemblée
  67. C HISTORIQUE : 09/04/04 : idmatp idmatd,
  68. C la cl dirichlet n'est plus forcément sur la
  69. C diagonale
  70. C HISTORIQUE :
  71. C***********************************************************************
  72. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  73. C en cas de modification de ce sous-programme afin de faciliter
  74. C la maintenance !
  75. C***********************************************************************
  76. -INC CCOPTIO
  77. -INC SMELEME
  78. POINTEUR ISPGT.MELEME
  79. -INC SMLENTI
  80. POINTEUR DDLPD.MLENTI
  81. POINTEUR DDLDP.MLENTI
  82. *-INC SMLLOGI
  83. * SEGMENT MLLOGI
  84. * LOGICAL LOGI(JG)
  85. * ENDSEGMENT
  86. * INTEGER JG
  87. * POINTEUR DDLCL.MLLOGI
  88. -INC SMCHPOI
  89. POINTEUR KCLIM.MCHPOI
  90. POINTEUR KMORS.PMORS
  91. POINTEUR KISA.IZA
  92. POINTEUR AMORS.PMORS
  93. POINTEUR AISA.IZA
  94. POINTEUR KS2B.IZA
  95. POINTEUR INCX.IZA
  96. POINTEUR IDMATP.IDMAT
  97. POINTEUR IDMATD.IDMAT
  98. *
  99. INTEGER IMPR,IRET
  100. *
  101. * .. Variables locales
  102. * .. Parameters
  103. REAL*8 ZERO,ONE
  104. PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  105. * ..
  106. CHARACTER*8 NOMINC
  107. LOGICAL FLINC
  108. LOGICAL FLDIA
  109. *
  110. INTEGER I1,IN,ISOUPO,INC,INBI,INBVA
  111. INTEGER N ,NSOUPO,NC ,NBI
  112. INTEGER ICOL,INTTT,JNTTT
  113. INTEGER NTTT
  114. REAL*8 TEMP
  115. LOGICAL LINCX
  116. C***
  117. LINCX=(INCX.NE.0)
  118. C On récupère les segments utiles
  119. IF (IMPR.GT.5) THEN
  120. WRITE(IOIMP,*) 'melim.eso : CL de Dirichlet'
  121. ENDIF
  122. SEGACT MATRIK
  123. MINC =KMINC
  124. ISPGT=KISPGT
  125. NTTT =KNTTT
  126. IDMATP=KIDMAT(1)
  127. IDMATD=KIDMAT(2)
  128. AMORS=KIDMAT(4)
  129. AISA =KIDMAT(5)
  130. SEGDES MATRIK
  131. C S'il n'y a pas de CL de Dirichlet, on n'a rien à faire
  132. IF(KCLIM.EQ.0) THEN
  133. SEGINI,KMORS=AMORS
  134. SEGDES KMORS
  135. SEGINI,KISA=AISA
  136. SEGDES KISA
  137. ELSE
  138. C 1) On s'occupe du second membre et de l'inconnue
  139. C On parcourt le chpoint un peu comme dans ch2vec.eso
  140. C Les DDL ou on applique une condition de Dirichlet
  141. C ont un DDLCL(No DDL)=1 et 0 sinon
  142. IF (IMPR.GT.6) THEN
  143. WRITE(IOIMP,*)
  144. $ 'Dirichlet sur l''inconnue et le second membre'
  145. ENDIF
  146. JG=NTTT
  147. SEGINI DDLPD
  148. SEGINI DDLDP
  149. SEGACT KS2B*MOD
  150. IF (LINCX) SEGACT INCX*MOD
  151. MCHPOI=KCLIM
  152. SEGACT MINC
  153. NBI=LISINC(/2)
  154. C In KRIPAD : SEGACT ISPGT
  155. C SEGINI MLENTI
  156. CALL KRIPAD(ISPGT,MLENTI)
  157. SEGDES ISPGT
  158. SEGACT IDMATP
  159. SEGACT IDMATD
  160. SEGACT MCHPOI
  161. NSOUPO=IPCHP(/1)
  162. DO 1 ISOUPO=1,NSOUPO
  163. MSOUPO=IPCHP(ISOUPO)
  164. SEGACT MSOUPO
  165. NC=NOCOMP(/2)
  166. MELEME=IGEOC
  167. MPOVAL=IPOVAL
  168. SEGACT MELEME
  169. SEGACT MPOVAL
  170. C NPTI=VPOCHA(/1)
  171. DO 11 INC=1,NC
  172. NOMINC=NOCOMP(INC)//' '
  173. FLINC=.FALSE.
  174. C Repeat..until
  175. INBI=1
  176. 111 CONTINUE
  177. IF (NOMINC.EQ.LISINC(INBI)) THEN
  178. FLINC=.TRUE.
  179. ELSEIF (INBI.LT.NBI) THEN
  180. INBI=INBI+1
  181. GOTO 111
  182. ENDIF
  183. IF (.NOT.FLINC) THEN
  184. * WRITE(IOIMP,*) 'melim.eso : composante ',NOMINC
  185. * $ ,'unknown'
  186. WRITE(IOIMP,*) 'CLIM : component name ',NOMINC
  187. $ ,' not in the matrix'
  188. GOTO 9999
  189. ELSE
  190. N=VPOCHA(/1)
  191. DO 112 IN=1,N
  192. I1=LECT(NUM(1,IN))
  193. IF(I1.EQ.0)THEN
  194. * WRITE(IOIMP,*) 'melim.eso : le point ',IN
  195. * $ ,' du ch.'
  196. * WRITE(IOIMP,*) 'n''appartient pas au vec.'
  197. WRITE(IOIMP,*) 'CLIM : point number ',(num(1,IN)),
  198. $ ' unknown ',NOMINC,' not in the matrix'
  199. GOTO 9999
  200. ELSE
  201. IF (MPOS(I1,INBI).EQ.0) THEN
  202. WRITE(IOIMP,*) 'CLIM : point number ',(num(1,IN)),
  203. $ ' unknown ',NOMINC,' not in the matrix'
  204. GOTO 9999
  205. ELSE
  206. IDIAGP=IDMATP.NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  207. IDIAGD=IDMATD.NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  208. C* INBVA=NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  209. C INPT=NUAN(I1)
  210. C INBVA=NPOS(INPT)+MPOS(INPT,INBI)-1
  211. C KVEC.A(INBVA)=KVEC.A(INBVA)+VPOCHA(IN,INC)
  212. C Je préfère surcharger
  213. C
  214. * KS2B.A(INBVA)=VPOCHA(IN,INC)
  215. * IF (LINCX) INCX.A(INBVA)=VPOCHA(IN,INC)
  216. KS2B.A(IDIAGD)=VPOCHA(IN,INC)
  217. IF (LINCX) INCX.A(IDIAGP)=VPOCHA(IN,INC)
  218. DDLPD.LECT(IDIAGP)=IDIAGD
  219. DDLDP.LECT(IDIAGD)=IDIAGP
  220. * DDLCL.LOGI(INBVA)=.TRUE.
  221. ENDIF
  222. ENDIF
  223. 112 CONTINUE
  224. ENDIF
  225. 11 CONTINUE
  226. SEGDES MSOUPO
  227. SEGDES MELEME
  228. SEGDES MPOVAL
  229. 1 CONTINUE
  230. SEGDES MCHPOI
  231. SEGDES IDMATD
  232. SEGDES IDMATP
  233. SEGSUP MLENTI
  234. SEGDES MINC
  235. IF (IMPR.GT.7) THEN
  236. SEGPRT,KS2B
  237. IF (LINCX) SEGPRT,INCX
  238. ENDIF
  239. IF (LINCX) SEGDES,INCX
  240. C 2) On s'occupe de la matrice :
  241. C pour les no de DDL ou on impose une CL dirichlet, il faut :
  242. C annuler la ligne correspondante sauf la diagonale qui vaut 1
  243. C reporter -colonne*clim sur le second membre
  244. C On parcourt donc toute la matrice
  245. IF (IMPR.GT.6) THEN
  246. WRITE(IOIMP,*) 'Dirichlet sur la matrice'
  247. ENDIF
  248. SEGINI,KMORS=AMORS
  249. SEGINI,KISA=AISA
  250. DO 2 INTTT=1,NTTT
  251. *! IF (DDLCL.LOGI(INTTT)) THEN
  252. IDIAGP=DDLDP.LECT(INTTT)
  253. IF (IDIAGP.NE.0) THEN
  254. C
  255. C CL de dirichlet pour la ligne INTTT :
  256. C On met 1 sur la diagonale et 0 ailleurs,
  257. C le second membre ne change pas
  258. C
  259. FLDIA=.FALSE.
  260. DO 21 ICOL=KMORS.IA(INTTT),(KMORS.IA(INTTT+1)-1)
  261. *! IF (KMORS.JA(ICOL).NE.INTTT) THEN
  262. IF (KMORS.JA(ICOL).NE.IDIAGP) THEN
  263. KISA.A(ICOL)=ZERO
  264. ELSE
  265. KISA.A(ICOL)=ONE
  266. FLDIA=.TRUE.
  267. ENDIF
  268. 21 CONTINUE
  269. IF (.NOT.FLDIA) THEN
  270. *! WRITE(IOIMP,*) 'melim.eso : diag.',INTTT,'inexistante'
  271. WRITE(IOIMP,*) 'melim.eso : diag.',IDIAGP,
  272. $ 'inexistante'
  273. GOTO 9999
  274. ENDIF
  275. ELSE
  276. C
  277. C Ligne INTTT sans condition de dirichlet :
  278. C On annule les termes des colonnes liées à une CL dirichlet,
  279. C et on reporte le produit colonne*CL au second membre
  280. C
  281. TEMP=ZERO
  282. DO 22 ICOL=KMORS.IA(INTTT),(KMORS.IA(INTTT+1)-1)
  283. JNTTT=KMORS.JA(ICOL)
  284. IDIAGD=DDLPD.LECT(JNTTT)
  285. *! IF (DDLCL.LOGI(JNTTT)) THEN
  286. IF (IDIAGD.NE.0) THEN
  287. *! IF (KS2B.A(JNTTT).NE.ZERO) THEN
  288. *! TEMP=TEMP+(KISA.A(ICOL)*KS2B.A(JNTTT))
  289. TEMP=TEMP+(KISA.A(ICOL)*KS2B.A(IDIAGD))
  290. *! ENDIF
  291. KISA.A(ICOL)=ZERO
  292. ENDIF
  293. 22 CONTINUE
  294. IF (TEMP.NE.ZERO) THEN
  295. KS2B.A(INTTT)=KS2B.A(INTTT)-TEMP
  296. ENDIF
  297. ENDIF
  298. 2 CONTINUE
  299. SEGDES KISA
  300. SEGDES KMORS
  301. SEGDES KS2B
  302. SEGSUP DDLPD
  303. SEGSUP DDLDP
  304. ENDIF
  305. *
  306. * Terminaison normale
  307. *
  308. IRET=0
  309. RETURN
  310. *
  311. * Format handling
  312. *
  313. 1002 FORMAT(10(1X,1PE11.4))
  314. *
  315. * Error handling
  316. *
  317. 9999 CONTINUE
  318. IRET=1
  319. WRITE(IOIMP,*) 'An error was detected in melim.eso'
  320. RETURN
  321. *
  322. * End of MELIM
  323. *
  324. END
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  
  341.  
  342.  
  343.  
  344.  
  345.  
  346.  

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