Télécharger matp.eso

Retour à la liste

Numérotation des lignes :

  1. C MATP SOURCE CB215821 19/08/20 21:19:35 10287
  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. CALL ACTOBJ(TYPOBJ,IOBRE,1)
  138. IF (IERR.NE.0) RETURN
  139. IPCK = IOBRE
  140. ELSE
  141. THETA = 0.D0
  142. DELTAT = 0.D0
  143. IPCK = 0
  144. ENDIF
  145. *
  146. * Test de l'ojet MODELE
  147. * Récupération des zones élémentaires ou DARCY est défini
  148. *
  149. SEGACT MMODEL
  150. NSOUS = KMODEL(/1)
  151. SEGINI IPMAHY
  152. IDARCY = 0
  153. IPT2=IPGEOM
  154. IPT1=IPGEOM
  155. DO 10 ISOUS=1,NSOUS
  156. IF(NSOUS.GT.1)THEN
  157. SEGACT IPT2
  158. IPT1=IPT2.LISOUS(ISOUS)
  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. 10 CONTINUE
  168. IF (IDARCY.EQ.0) THEN
  169. MOTERR = LETYPE
  170. CALL ERREUR(193)
  171. GOTO 100
  172. ENDIF
  173. *
  174. * Recuperation des pointeurs ELTFA pour les zones ou DARCY est defini
  175. *
  176. MELEME = IELTFA
  177. SEGACT MELEME
  178. LZONES = LISOUS(/1)
  179. IF (LZONES.EQ.0) LZONES = 1
  180. IPT1 = IPGEOM
  181. SEGACT IPT1
  182. DO 30 ISOUS=1,NSOUS
  183. IMAMEL = MAHYBR(ISOUS)
  184. IF (IMAMEL.NE.0) THEN
  185. DO 20 ISZ=1,LZONES
  186. IF (LZONES.EQ.1) THEN
  187. IPT2 = IPT1
  188. IPT3 = MELEME
  189. ELSE
  190. IPT2 = IPT1.LISOUS(ISZ)
  191. IPT3 = LISOUS(ISZ)
  192. ENDIF
  193. IF (IPT2.EQ.IMAMEL) THEN
  194. MAHYBR(ISOUS) = IPT3
  195. GOTO 30
  196. ENDIF
  197. 20 CONTINUE
  198. IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN
  199. MOTERR(1:8) = ' MODELE '
  200. MOTERR(9:16) = ' ELTFA '
  201. INTERR(1) = ISOUS
  202. CALL ERREUR(698)
  203. GOTO 100
  204. ENDIF
  205. ENDIF
  206. 30 CONTINUE
  207. *
  208. * Test du sous-type de la matrice de rigidité récupérée
  209. * Controle des MELEME de l'objet RIGIDITE % au ELTFA liée au MMODEL
  210. *
  211. SEGACT MRIGID
  212. LETYPE = MTYMAT
  213. IF (LETYPE.NE.'DARCY') THEN
  214. MOTERR(1:8) = 'RIGIDITE'
  215. MOTERR(9:16) = 'DARCY'
  216. CALL ERREUR(79)
  217. SEGDES MRIGID
  218. GOTO 100
  219. ENDIF
  220. DO 40 ISOUS=1,NSOUS
  221. IF (MAHYBR(ISOUS).NE.0) THEN
  222. IPTEST = IRIGEL(1,ISOUS)
  223. IF (IPTEST.NE.MAHYBR(ISOUS)) THEN
  224. MOTERR(1:8) = 'DARCY'
  225. MOTERR(9:16) = ' ELTFA '
  226. INTERR(1) = ISOUS
  227. CALL ERREUR(698)
  228. SEGDES MRIGID
  229. GOTO 100
  230. ENDIF
  231. ENDIF
  232. 40 CONTINUE
  233. SEGDES MRIGID
  234. *
  235. * Vérification du support du MCHAML % au MMODEL
  236. * Controle du nombre de composantes reelles
  237. *
  238. IF (IPCK.NE.0) THEN
  239. ISUP = 2
  240. ICOND = 1
  241. CALL QUESUP(IPMODE,IPCK,ISUP,ICOND,IRET,IRET2)
  242. IF (IRET.GT.1) GOTO 100
  243. MCHELM = IPCK
  244. SEGACT MCHELM
  245. DO 50 ISOUS=1,NSOUS
  246. IF (MAHYBR(ISOUS).NE.0) THEN
  247. MCHAML = ICHAML(ISOUS)
  248. SEGACT MCHAML
  249. N2 = TYPCHE(/2)
  250. IF (N2.GT.1) THEN
  251. CALL ERREUR(320)
  252. GOTO 100
  253. ENDIF
  254. CHAR6 = TYPCHE(1)(1:6)
  255. IF (CHAR6.NE.'REAL*8') THEN
  256. MOTERR(1:8) = NOMCHE(1)
  257. CALL ERREUR(679)
  258. GOTO 100
  259. ENDIF
  260. ENDIF
  261. 50 CONTINUE
  262. ENDIF
  263. SEGDES IPMAHY
  264. *
  265. * Construction de IPRIG1
  266. *
  267. CALL MATP1(NSOUS,IPGEOM,IPMAHY,IPRIGI,THETA,DELTAT,IPCK,IPRIG1)
  268. CALL ECROBJ('RIGIDITE',IPRIG1)
  269. *
  270. * Ménage
  271. *
  272. 100 CONTINUE
  273. SEGSUP IPMAHY
  274. END
  275.  
  276.  
  277.  

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