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.  
  77. -INC PPARAM
  78. -INC CCOPTIO
  79. -INC SMELEME
  80. POINTEUR ISPGT.MELEME
  81. -INC SMLENTI
  82. POINTEUR DDLPD.MLENTI
  83. POINTEUR DDLDP.MLENTI
  84. *-INC SMLLOGI
  85. * SEGMENT MLLOGI
  86. * LOGICAL LOGI(JG)
  87. * ENDSEGMENT
  88. * INTEGER JG
  89. * POINTEUR DDLCL.MLLOGI
  90. -INC SMCHPOI
  91. POINTEUR KCLIM.MCHPOI
  92. POINTEUR KMORS.PMORS
  93. POINTEUR KISA.IZA
  94. POINTEUR AMORS.PMORS
  95. POINTEUR AISA.IZA
  96. POINTEUR KS2B.IZA
  97. POINTEUR INCX.IZA
  98. POINTEUR IDMATP.IDMAT
  99. POINTEUR IDMATD.IDMAT
  100. *
  101. INTEGER IMPR,IRET
  102. *
  103. * .. Variables locales
  104. * .. Parameters
  105. REAL*8 ZERO,ONE
  106. PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  107. * ..
  108. CHARACTER*8 NOMINC
  109. LOGICAL FLINC
  110. LOGICAL FLDIA
  111. *
  112. INTEGER I1,IN,ISOUPO,INC,INBI,INBVA
  113. INTEGER N ,NSOUPO,NC ,NBI
  114. INTEGER ICOL,INTTT,JNTTT
  115. INTEGER NTTT
  116. REAL*8 TEMP
  117. LOGICAL LINCX
  118. C***
  119. LINCX=(INCX.NE.0)
  120. C On récupère les segments utiles
  121. IF (IMPR.GT.5) THEN
  122. WRITE(IOIMP,*) 'melim.eso : CL de Dirichlet'
  123. ENDIF
  124. SEGACT MATRIK
  125. MINC =KMINC
  126. ISPGT=KISPGT
  127. NTTT =KNTTT
  128. IDMATP=KIDMAT(1)
  129. IDMATD=KIDMAT(2)
  130. AMORS=KIDMAT(4)
  131. AISA =KIDMAT(5)
  132. SEGDES MATRIK
  133. C S'il n'y a pas de CL de Dirichlet, on n'a rien à faire
  134. IF(KCLIM.EQ.0) THEN
  135. SEGINI,KMORS=AMORS
  136. SEGDES KMORS
  137. SEGINI,KISA=AISA
  138. SEGDES KISA
  139. ELSE
  140. C 1) On s'occupe du second membre et de l'inconnue
  141. C On parcourt le chpoint un peu comme dans ch2vec.eso
  142. C Les DDL ou on applique une condition de Dirichlet
  143. C ont un DDLCL(No DDL)=1 et 0 sinon
  144. IF (IMPR.GT.6) THEN
  145. WRITE(IOIMP,*)
  146. $ 'Dirichlet sur l''inconnue et le second membre'
  147. ENDIF
  148. JG=NTTT
  149. SEGINI DDLPD
  150. SEGINI DDLDP
  151. SEGACT KS2B*MOD
  152. IF (LINCX) SEGACT INCX*MOD
  153. MCHPOI=KCLIM
  154. SEGACT MINC
  155. NBI=LISINC(/2)
  156. C In KRIPAD : SEGACT ISPGT
  157. C SEGINI MLENTI
  158. CALL KRIPAD(ISPGT,MLENTI)
  159. SEGDES ISPGT
  160. SEGACT IDMATP
  161. SEGACT IDMATD
  162. SEGACT MCHPOI
  163. NSOUPO=IPCHP(/1)
  164. DO 1 ISOUPO=1,NSOUPO
  165. MSOUPO=IPCHP(ISOUPO)
  166. SEGACT MSOUPO
  167. NC=NOCOMP(/2)
  168. MELEME=IGEOC
  169. MPOVAL=IPOVAL
  170. SEGACT MELEME
  171. SEGACT MPOVAL
  172. C NPTI=VPOCHA(/1)
  173. DO 11 INC=1,NC
  174. NOMINC=NOCOMP(INC)//' '
  175. FLINC=.FALSE.
  176. C Repeat..until
  177. INBI=1
  178. 111 CONTINUE
  179. IF (NOMINC.EQ.LISINC(INBI)) THEN
  180. FLINC=.TRUE.
  181. ELSEIF (INBI.LT.NBI) THEN
  182. INBI=INBI+1
  183. GOTO 111
  184. ENDIF
  185. IF (.NOT.FLINC) THEN
  186. * WRITE(IOIMP,*) 'melim.eso : composante ',NOMINC
  187. * $ ,'unknown'
  188. WRITE(IOIMP,*) 'CLIM : component name ',NOMINC
  189. $ ,' not in the matrix'
  190. GOTO 9999
  191. ELSE
  192. N=VPOCHA(/1)
  193. DO 112 IN=1,N
  194. I1=LECT(NUM(1,IN))
  195. IF(I1.EQ.0)THEN
  196. * WRITE(IOIMP,*) 'melim.eso : le point ',IN
  197. * $ ,' du ch.'
  198. * WRITE(IOIMP,*) 'n''appartient pas au vec.'
  199. WRITE(IOIMP,*) 'CLIM : point number ',(num(1,IN)),
  200. $ ' unknown ',NOMINC,' not in the matrix'
  201. GOTO 9999
  202. ELSE
  203. IF (MPOS(I1,INBI).EQ.0) THEN
  204. WRITE(IOIMP,*) 'CLIM : point number ',(num(1,IN)),
  205. $ ' unknown ',NOMINC,' not in the matrix'
  206. GOTO 9999
  207. ELSE
  208. IDIAGP=IDMATP.NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  209. IDIAGD=IDMATD.NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  210. C* INBVA=NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  211. C INPT=NUAN(I1)
  212. C INBVA=NPOS(INPT)+MPOS(INPT,INBI)-1
  213. C KVEC.A(INBVA)=KVEC.A(INBVA)+VPOCHA(IN,INC)
  214. C Je préfère surcharger
  215. C
  216. * KS2B.A(INBVA)=VPOCHA(IN,INC)
  217. * IF (LINCX) INCX.A(INBVA)=VPOCHA(IN,INC)
  218. KS2B.A(IDIAGD)=VPOCHA(IN,INC)
  219. IF (LINCX) INCX.A(IDIAGP)=VPOCHA(IN,INC)
  220. DDLPD.LECT(IDIAGP)=IDIAGD
  221. DDLDP.LECT(IDIAGD)=IDIAGP
  222. * DDLCL.LOGI(INBVA)=.TRUE.
  223. ENDIF
  224. ENDIF
  225. 112 CONTINUE
  226. ENDIF
  227. 11 CONTINUE
  228. SEGDES MSOUPO
  229. SEGDES MELEME
  230. SEGDES MPOVAL
  231. 1 CONTINUE
  232. SEGDES MCHPOI
  233. SEGDES IDMATD
  234. SEGDES IDMATP
  235. SEGSUP MLENTI
  236. SEGDES MINC
  237. IF (IMPR.GT.7) THEN
  238. SEGPRT,KS2B
  239. IF (LINCX) SEGPRT,INCX
  240. ENDIF
  241. IF (LINCX) SEGDES,INCX
  242. C 2) On s'occupe de la matrice :
  243. C pour les no de DDL ou on impose une CL dirichlet, il faut :
  244. C annuler la ligne correspondante sauf la diagonale qui vaut 1
  245. C reporter -colonne*clim sur le second membre
  246. C On parcourt donc toute la matrice
  247. IF (IMPR.GT.6) THEN
  248. WRITE(IOIMP,*) 'Dirichlet sur la matrice'
  249. ENDIF
  250. SEGINI,KMORS=AMORS
  251. SEGINI,KISA=AISA
  252. DO 2 INTTT=1,NTTT
  253. *! IF (DDLCL.LOGI(INTTT)) THEN
  254. IDIAGP=DDLDP.LECT(INTTT)
  255. IF (IDIAGP.NE.0) THEN
  256. C
  257. C CL de dirichlet pour la ligne INTTT :
  258. C On met 1 sur la diagonale et 0 ailleurs,
  259. C le second membre ne change pas
  260. C
  261. FLDIA=.FALSE.
  262. DO 21 ICOL=KMORS.IA(INTTT),(KMORS.IA(INTTT+1)-1)
  263. *! IF (KMORS.JA(ICOL).NE.INTTT) THEN
  264. IF (KMORS.JA(ICOL).NE.IDIAGP) THEN
  265. KISA.A(ICOL)=ZERO
  266. ELSE
  267. KISA.A(ICOL)=ONE
  268. FLDIA=.TRUE.
  269. ENDIF
  270. 21 CONTINUE
  271. IF (.NOT.FLDIA) THEN
  272. *! WRITE(IOIMP,*) 'melim.eso : diag.',INTTT,'inexistante'
  273. WRITE(IOIMP,*) 'melim.eso : diag.',IDIAGP,
  274. $ 'inexistante'
  275. GOTO 9999
  276. ENDIF
  277. ELSE
  278. C
  279. C Ligne INTTT sans condition de dirichlet :
  280. C On annule les termes des colonnes liées à une CL dirichlet,
  281. C et on reporte le produit colonne*CL au second membre
  282. C
  283. TEMP=ZERO
  284. DO 22 ICOL=KMORS.IA(INTTT),(KMORS.IA(INTTT+1)-1)
  285. JNTTT=KMORS.JA(ICOL)
  286. IDIAGD=DDLPD.LECT(JNTTT)
  287. *! IF (DDLCL.LOGI(JNTTT)) THEN
  288. IF (IDIAGD.NE.0) THEN
  289. *! IF (KS2B.A(JNTTT).NE.ZERO) THEN
  290. *! TEMP=TEMP+(KISA.A(ICOL)*KS2B.A(JNTTT))
  291. TEMP=TEMP+(KISA.A(ICOL)*KS2B.A(IDIAGD))
  292. *! ENDIF
  293. KISA.A(ICOL)=ZERO
  294. ENDIF
  295. 22 CONTINUE
  296. IF (TEMP.NE.ZERO) THEN
  297. KS2B.A(INTTT)=KS2B.A(INTTT)-TEMP
  298. ENDIF
  299. ENDIF
  300. 2 CONTINUE
  301. SEGDES KISA
  302. SEGDES KMORS
  303. SEGDES KS2B
  304. SEGSUP DDLPD
  305. SEGSUP DDLDP
  306. ENDIF
  307. *
  308. * Terminaison normale
  309. *
  310. IRET=0
  311. RETURN
  312. *
  313. * Format handling
  314. *
  315. 1002 FORMAT(10(1X,1PE11.4))
  316. *
  317. * Error handling
  318. *
  319. 9999 CONTINUE
  320. IRET=1
  321. WRITE(IOIMP,*) 'An error was detected in melim.eso'
  322. RETURN
  323. *
  324. * End of MELIM
  325. *
  326. END
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  
  341.  
  342.  
  343.  
  344.  
  345.  
  346.  
  347.  
  348.  

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