Télécharger matp.eso

Retour à la liste

Numérotation des lignes :

  1. C MATP SOURCE CHAT 11/03/16 21:27:48 6902
  2. SUBROUTINE MATP
  3. C-----------------------------------------------------------------------
  4. C Cette subroutine permet de creer les matrices elementaires du système
  5. C matriciel en trace de charge dans le cadre de la résolution des
  6. C équations de DARCY par une méthode d'éléments finis mixtes hybride.
  7. C On obtient ce systeme en exprimant les inconnues de charge et de
  8. C flux en fonction des traces de charge.
  9. C -1 t -1
  10. C M1 = M2 * ( COEF * DIV * DIV * M2 - Id )
  11. C
  12. C avec M2 : Matrice massse hybride
  13. C DIV : Matrice uniligne representant la divergence
  14. C COEF : Scalaire valant 1/d dans le cas permanant, THETA*BETA/d
  15. C dans le cas transitoire.
  16. C -1 t
  17. C d : Scalaire valant ( DIV * M2 * DIV )
  18. C
  19. C THETA : Scalaire de la THETA-méthode | Données
  20. C BETA : Scalaire d*dt/(Ck|K| + THETA*d*dt) | utiles
  21. C dt : Scalaire pas de temps | uniquement
  22. C Ck : Scalaire coefficient d'emmagasinement | en
  23. C |K| : Scalaire aire de l'élément K | transitoire
  24. C
  25. C---------------------------
  26. C Phrase d'appel (GIBIANE) :
  27. C---------------------------
  28. C
  29. C MRIGI1 = MATP MMODEL MRIGI2 (TABLE2) ;
  30. C
  31. C------------------------
  32. C Operandes et resultat :
  33. C------------------------
  34. C
  35. C MRIGI1 : Matrices elementaires en trace de charge.
  36. C MMODEL : Objet modele specifiant la formulation.
  37. C MRIGI2 : Matrices masses hybrides elementaires creees par MHYB.
  38. C TABLE2 : TABLE DARCY_TRANSITOIRE contenant les infos transitoires.
  39. C
  40. C-----------------------------------------------------------------------
  41. C
  42. C Langage : ESOPE + FORTRAN77
  43. C
  44. C Auteurs : 09/93 F.DABBENE - Cas permanent
  45. C 09/94 X.NOUVELLON - Extension au cas transitoire
  46. C
  47. C-----------------------------------------------------------------------
  48. IMPLICIT INTEGER(I-N)
  49. IMPLICIT REAL*8 (A-H,O-Z)
  50. *
  51. -INC CCOPTIO
  52. -INC SMELEME
  53. -INC SMRIGID
  54. -INC SMMODEL
  55. -INC SMTABLE
  56. -INC SMCHAML
  57. *
  58. SEGMENT IPMAHY
  59. INTEGER MAHYBR(NSOUS)
  60. ENDSEGMENT
  61. *
  62. LOGICAL LOGRE,LOGIN
  63. CHARACTER*6 CHAR6
  64. CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE
  65. *
  66. * Initialisations
  67. *
  68. IVALIN = 0
  69. XVALIN = 0.D0
  70. LOGIN = .TRUE.
  71. IOBIN = 0
  72. TAPIND = 'MOT '
  73. *
  74. * Lecture du MMODEL
  75. *
  76. CALL LIROBJ('MMODEL',IPMODE,1,IRET1)
  77. IF (IERR.NE.0) RETURN
  78. MMODEL = IPMODE
  79. *
  80. * Lecture de la TABLE domaine et récupération des infos
  81. *
  82. IPTABL = 0
  83. CALL LEKMOD(IPMODE,IPTABL,IRET)
  84. IF (IERR.NE.0) RETURN
  85. CHARIN = 'MAILLAGE'
  86. TYPOBJ = 'MAILLAGE'
  87. CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  88. IF (IERR.NE.0) RETURN
  89. IPGEOM = IOBRE
  90. CALL LEKTAB(IPTABL,'ELTFA',IOBRE)
  91. IF (IERR.NE.0) RETURN
  92. IELTFA = IOBRE
  93. *
  94. * Lecture de RIGIDITE
  95. *
  96. CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET1)
  97. IF (IERR.NE.0) RETURN
  98. MRIGID = IPRIGI
  99. *
  100. * Lecture éventuelle de la TABLE DARCY_TRANSITOIRE
  101. *
  102. IPTABL = 0
  103. CALL LIRTAB('DARCY_TRANSITOIRE',IPTABL,0,IRET)
  104. IF (IERR.NE.0) RETURN
  105. IF (IRET.EQ.1) THEN
  106. TYPOBJ = 'FLOTTANT'
  107. CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'THETA',LOGIN,IOBIN,
  108. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  109. IF (IERR.NE.0) RETURN
  110. THETA = XVALRE
  111. THEMIN = -1.D-12
  112. THEMAX = 1.D0 - THEMIN
  113. IF ((THETA.LT.THEMIN).OR.(THETA.GT.THEMAX)) THEN
  114. REAERR(1) = REAL(THETA)
  115. REAERR(2) = REAL(0.D0)
  116. REAERR(3) = REAL(1.D0)
  117. MOTERR(1:8) = ' THETA '
  118. CALL ERREUR(42)
  119. RETURN
  120. ENDIF
  121. IF (THETA.LT.0.D0) THETA=0.D0
  122. IF (THETA.GT.1.D0) THETA=1.D0
  123. CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'PAS',LOGIN,IOBIN,
  124. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  125. IF (IERR.NE.0) RETURN
  126. DELTAT = XVALRE
  127. IF (DELTAT.LT.0.D0) THEN
  128. REAERR(1) = REAL(DELTAT)
  129. REAERR(2) = REAL(0.D0)
  130. MOTERR(1:8) = ' DELTAT '
  131. CALL ERREUR(41)
  132. RETURN
  133. ENDIF
  134. TYPOBJ = 'MCHAML '
  135. CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'SURF',LOGIN,IOBIN,
  136. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  137. IF (IERR.NE.0) RETURN
  138. IPCK = IOBRE
  139. ELSE
  140. THETA = 0.D0
  141. DELTAT = 0.D0
  142. IPCK = 0
  143. ENDIF
  144. *
  145. * Test de l'ojet MODELE
  146. * Récupération des zones élémentaires ou DARCY est défini
  147. *
  148. SEGACT MMODEL
  149. NSOUS = KMODEL(/1)
  150. SEGINI IPMAHY
  151. IDARCY = 0
  152. IPT2=IPGEOM
  153. IPT1=IPGEOM
  154. DO 10 ISOUS=1,NSOUS
  155. IF(NSOUS.GT.1)THEN
  156. SEGACT IPT2
  157. IPT1=IPT2.LISOUS(ISOUS)
  158. SEGDES IPT2
  159. ENDIF
  160. IMODEL = KMODEL(ISOUS)
  161. SEGACT IMODEL
  162. LETYPE = FORMOD(1)
  163. IF (LETYPE.EQ.'DARCY') THEN
  164. IDARCY = IDARCY + 1
  165. MAHYBR(ISOUS) = IPT1
  166. ENDIF
  167. SEGDES IMODEL
  168. 10 CONTINUE
  169. SEGDES MMODEL
  170. IF (IDARCY.EQ.0) THEN
  171. MOTERR = LETYPE
  172. CALL ERREUR(193)
  173. GOTO 100
  174. ENDIF
  175. *
  176. * Recuperation des pointeurs ELTFA pour les zones ou DARCY est defini
  177. *
  178. MELEME = IELTFA
  179. SEGACT MELEME
  180. LZONES = LISOUS(/1)
  181. IF (LZONES.EQ.0) LZONES = 1
  182. IPT1 = IPGEOM
  183. SEGACT IPT1
  184. DO 30 ISOUS=1,NSOUS
  185. IMAMEL = MAHYBR(ISOUS)
  186. IF (IMAMEL.NE.0) THEN
  187. DO 20 ISZ=1,LZONES
  188. IF (LZONES.EQ.1) THEN
  189. IPT2 = IPT1
  190. IPT3 = MELEME
  191. ELSE
  192. IPT2 = IPT1.LISOUS(ISZ)
  193. IPT3 = LISOUS(ISZ)
  194. ENDIF
  195. IF (IPT2.EQ.IMAMEL) THEN
  196. MAHYBR(ISOUS) = IPT3
  197. GOTO 30
  198. ENDIF
  199. 20 CONTINUE
  200. IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN
  201. MOTERR(1:8) = ' MODELE '
  202. MOTERR(9:16) = ' ELTFA '
  203. INTERR(1) = ISOUS
  204. CALL ERREUR(698)
  205. SEGDES IPT1
  206. SEGDES MELEME
  207. GOTO 100
  208. ENDIF
  209. ENDIF
  210. 30 CONTINUE
  211. SEGDES IPT1
  212. SEGDES MELEME
  213. *
  214. * Test du sous-type de la matrice de rigidité récupérée
  215. * Controle des MELEME de l'objet RIGIDITE % au ELTFA liée au MMODEL
  216. *
  217. SEGACT MRIGID
  218. LETYPE = MTYMAT
  219. IF (LETYPE.NE.'DARCY') THEN
  220. MOTERR(1:8) = 'RIGIDITE'
  221. MOTERR(9:16) = 'DARCY'
  222. CALL ERREUR(79)
  223. SEGDES MRIGID
  224. GOTO 100
  225. ENDIF
  226. DO 40 ISOUS=1,NSOUS
  227. IF (MAHYBR(ISOUS).NE.0) THEN
  228. IPTEST = IRIGEL(1,ISOUS)
  229. IF (IPTEST.NE.MAHYBR(ISOUS)) THEN
  230. MOTERR(1:8) = 'DARCY'
  231. MOTERR(9:16) = ' ELTFA '
  232. INTERR(1) = ISOUS
  233. CALL ERREUR(698)
  234. SEGDES MRIGID
  235. GOTO 100
  236. ENDIF
  237. ENDIF
  238. 40 CONTINUE
  239. SEGDES MRIGID
  240. *
  241. * Vérification du support du MCHAML % au MMODEL
  242. * Controle du nombre de composantes reelles
  243. *
  244. IF (IPCK.NE.0) THEN
  245. ISUP = 2
  246. ICOND = 1
  247. CALL QUESUP(IPMODE,IPCK,ISUP,ICOND,IRET,IRET2)
  248. IF (IRET.GT.1) GOTO 100
  249. MCHELM = IPCK
  250. SEGACT MCHELM
  251. DO 50 ISOUS=1,NSOUS
  252. IF (MAHYBR(ISOUS).NE.0) THEN
  253. MCHAML = ICHAML(ISOUS)
  254. SEGACT MCHAML
  255. N2 = TYPCHE(/2)
  256. IF (N2.GT.1) THEN
  257. CALL ERREUR(320)
  258. SEGDES MCHELM
  259. SEGDES MCHAML
  260. GOTO 100
  261. ENDIF
  262. CHAR6 = TYPCHE(1)(1:6)
  263. IF (CHAR6.NE.'REAL*8') THEN
  264. MOTERR(1:8) = NOMCHE(1)
  265. CALL ERREUR(679)
  266. SEGDES MCHELM
  267. SEGDES MCHAML
  268. GOTO 100
  269. ENDIF
  270. SEGDES MCHAML
  271. ENDIF
  272. 50 CONTINUE
  273. SEGDES MCHELM
  274. ENDIF
  275. SEGDES IPMAHY
  276. *
  277. * Construction de IPRIG1
  278. *
  279. CALL MATP1(NSOUS,IPGEOM,IPMAHY,IPRIGI,THETA,DELTAT,IPCK,IPRIG1)
  280. CALL ECROBJ('RIGIDITE',IPRIG1)
  281. *
  282. * Ménage
  283. *
  284. 100 CONTINUE
  285. SEGSUP IPMAHY
  286. RETURN
  287. END
  288.  
  289.  
  290.  
  291.  
  292.  
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299.  
  300.  
  301.  
  302.  

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