Télécharger cadgsi.eso

Retour à la liste

Numérotation des lignes :

  1. C CADGSI SOURCE MAGN 17/02/24 21:15:02 9323
  2. SUBROUTINE CADGSI
  3. C************************************************************************
  4. C
  5. C OBJET :
  6. C
  7. C CALCUL DE LA MATRICE MASSE DIAGONALE ---> Creation d'un CHAMPOIN
  8. C Pour les SOMMETS D0=NI ( MASSE LUMPE )
  9. C Pour les FACES D0=1/2(Vol1 + Vol2)
  10. C Pour les CENTRES D0=Vol Elt
  11. C
  12. C SYNTAXE :
  13. C
  14. C RES = DGSI OBJ1 <TYPE> <'IMPR'> ;
  15. C
  16. C OBJ1 : Table DOMAINE
  17. C TYPE ; SOMMET , FACE , CENTRE (par defaut SOMMET) MSOMMET CENTREP1
  18. C
  19. C
  20. C************************************************************************
  21. IMPLICIT INTEGER(I-N)
  22. IMPLICIT REAL*8 (A-H,O-Z)
  23. -INC CCGEOME
  24. -INC SMCHAML
  25. -INC CCOPTIO
  26. -INC SMELEME
  27. POINTEUR MELEMP.MELEME
  28. -INC SMTABLE
  29. POINTEUR MTABD.MTABLE
  30. -INC SMCOORD
  31. -INC SMLENTI
  32. -INC SMCHPOI
  33. -INC SIZFFB
  34. POINTEUR IZF1.IZFFM,IZH2.IZHR
  35. CHARACTER*8 NOM0,CHAI,LISMO(5),TYPE,TYPC,NOM
  36. DATA LISMO/'SOMMET ','FACE ','CENTRE ','MSOMMET ','CENTREP1'/
  37. C ***************************************************************
  38. C----------------------------------------------------------------------
  39. C KPOIN = 0->SOMMET 1-> FACE 2-> CENTRE 3-> CENTREP0 4-> CENTREP1 5-> MSOMMET
  40. C INEFMD : Type formulation INEFMD=1 LINE,=2 MACRO,=3 QUADRATIQUE, INEFMD=4 LINB
  41. C************************************************************************
  42.  
  43.  
  44. IKAS=1
  45. IMPR=0
  46. IAXI=0
  47. IF(IFOMOD.EQ.0)IAXI=2
  48.  
  49. CALL LITABS('DOMAINE ',MTABD,1,1,IRET)
  50. IF(IRET.EQ.0)THEN
  51. WRITE(6,*)' On attend une table de soustype DOMAINE'
  52. RETURN
  53. ENDIF
  54.  
  55. 19 CONTINUE
  56. CALL LIRCHA(CHAI,0,IRET)
  57. IF(IRET.EQ.0)GO TO 20
  58. CALL OPTLI(IP,LISMO,CHAI,5)
  59. IF(IP.EQ.0)THEN
  60. WRITE(6,*)' On attend un mot cle parmi SOMMET FACE CENTRE ',
  61. & 'MSOMMET CENTREP1 '
  62. RETURN
  63. ENDIF
  64. IKAS=IP
  65.  
  66. 20 CONTINUE
  67.  
  68. IF(IKAS.EQ.1.OR.IKAS.EQ.4)THEN
  69. C SOMMET et MSOMMET
  70. C SOMMET
  71. IF(IKAS.EQ.1)THEN
  72. TYPE=' '
  73. CALL ACMO(MTABD,'MAILLAGE',TYPE,MELEME)
  74. IF(TYPE.NE.'MAILLAGE')RETURN
  75. TYPE=' '
  76. CALL ACMO(MTABD,'SOMMET',TYPE,MELEMS)
  77. IF(TYPE.NE.'MAILLAGE')RETURN
  78. CALL CRCHPT('SOMMET',MELEMS,1,MCHPOI)
  79. ENDIF
  80. C MSOMMET
  81. IF(IKAS.EQ.4)THEN
  82. c write(6,*)' CADGSI: IKAS=4'
  83. TYPE=' '
  84. CALL ACMO(MTABD,'MMAIL ',TYPE,MELEME)
  85. IF(TYPE.NE.'MAILLAGE')RETURN
  86. TYPE=' '
  87. CALL ACMO(MTABD,'MSOMMET',TYPE,MELEMS)
  88. IF(TYPE.NE.'MAILLAGE')RETURN
  89. CALL CRCHPT('MSOMMET',MELEMS,1,MCHPOI)
  90. ENDIF
  91.  
  92. CALL KRIPAD(MELEMS,MLENTI)
  93.  
  94. C CREATION DE LA DIAGONALE
  95. CALL LICHTM(MCHPOI,MPOVAL,TYPC,IGEOM)
  96. SEGACT MELEME
  97. NBSOUS=LISOUS(/1)
  98. IF(NBSOUS.EQ.0)NBSOUS=1
  99. C
  100. C BOUCLE SUR LES TYPES D'ELEMENTS ET CALCUL
  101. C
  102. DO 1 KSOUS=1,NBSOUS
  103. IF(NBSOUS.EQ.1)IPT1=MELEME
  104. IF(NBSOUS.GT.1)IPT1=LISOUS(KSOUS)
  105. SEGACT IPT1
  106.  
  107. NP=IPT1.NUM(/1)
  108. NEL=IPT1.NUM(/2)
  109. C
  110. NOM0=NOMS(IPT1.ITYPEL)//' '
  111. IF(INEFMD.EQ.1.AND.IKAS.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'P1P1'
  112. IF(INEFMD.EQ.2.AND.IKAS.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'MCF1'
  113. IF(INEFMD.EQ.3.AND.IKAS.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'PFP1'
  114. c write(6,*)'4 NOM0=',NOM0
  115. CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  116. SEGACT IZFFM*MOD
  117. IZHR=KZHR(1)
  118. SEGACT IZHR*MOD
  119. C
  120. NPG=FN(/2)
  121. NES=GR(/1)
  122.  
  123. IF(IMPR.NE.0)THEN
  124. WRITE(6,*)' SUB CADGSI : NES,NP,NPG,IDIM,NEL='
  125. & ,NES,NP,NPG,IDIM,NEL
  126. ENDIF
  127. C
  128. DO 10 K=1,NEL
  129. C
  130. NPGR=0
  131. IF(IAXI.NE.0)NPGR=NPG
  132. C
  133. DO 12 I=1,NP
  134. J=IPT1.NUM(I,K)
  135. DO 12 N=1,IDIM
  136. XYZ(N,I)=XCOOR((J-1)*(IDIM+1) +N)
  137. 12 CONTINUE
  138.  
  139. CALL CALJBC(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE)
  140.  
  141. IF(IMPR.NE.0)THEN
  142. WRITE(6,*)' SUB CADGSI : AIRE=',AIRE
  143. WRITE(6,*)' SUB CADGSI : LER '
  144. WRITE(6,1001)(IPT1.NUM(I,K),I=1,NP)
  145. WRITE(6,*)' SUB CADGSI : XYZ '
  146. WRITE(6,1002)((XYZ(N,I),N=1,2),I=1,NP)
  147. ENDIF
  148. C
  149. DO 3 J=1,NP
  150. SJ=0.D0
  151. DO 4 L=1,NPG
  152. SJ=SJ+FN(J,L)*PGSQ(L)
  153. 4 CONTINUE
  154. C SD(J,K)=SJ
  155. JU=LECT(IPT1.NUM(J,K))
  156. C D0(JU)=D0(JU)+SJ
  157. VPOCHA(JU,1)=VPOCHA(JU,1)+SJ
  158. 3 CONTINUE
  159. C
  160. 10 CONTINUE
  161.  
  162. SEGDES IPT1
  163. SEGSUP IZFFM,IZHR
  164. 1 CONTINUE
  165. SEGDES MELEME
  166. SEGSUP MLENTI
  167. SEGDES MPOVAL
  168.  
  169. IF(IMPR.NE.0)THEN
  170. WRITE(6,*)' SUB CADGSI : CALCUL DE LA DIAGONALE'
  171. C WRITE(6,1003)(I,VPOCHA(I,1),I=1,NPT)
  172. WRITE(6,*)' FIN DE CADGSI'
  173. ENDIF
  174.  
  175. C FACE
  176. ELSEIF(IKAS.EQ.2)THEN
  177.  
  178. TYPE=' '
  179. CALL ACMO(MTABD,'FACE',TYPE,MELEMF)
  180. IF(TYPE.NE.'MAILLAGE')RETURN
  181. TYPE=' '
  182. CALL ACMO(MTABD,'CENTRE',TYPE,MELEMC)
  183. IF(TYPE.NE.'MAILLAGE')RETURN
  184. CALL LEKTAB(MTABD,'FACEL',MELEME)
  185. IF(MELEME.EQ.0)RETURN
  186. CALL LEKTAB(MTABD,'XXVOLUM',MCHPO1)
  187. IF(MCHPO1.EQ.0)RETURN
  188.  
  189. CALL KRIPAD(MELEMC,MLENT1)
  190. CALL KRIPAD(MELEMF,MLENT2)
  191.  
  192. C CREATION DE LA DIAGONALE
  193. CALL CRCHPT('FACE',MELEMF,1,MCHPOI)
  194. CALL LICHTM(MCHPOI,MPOVAL,TYPC,IGEOM)
  195. CALL LICHTM(MCHPO1,MPOVA1,TYPC,IGEOM)
  196. SEGACT MELEME
  197. C
  198. C BOUCLE SUR LES TYPES D'ELEMENTS ET CALCUL
  199. C
  200. NBEL=NUM(/2)
  201. DO 2 K=1,NBEL
  202. I1=NUM(1,K)
  203. I2=NUM(2,K)
  204. I3=NUM(3,K)
  205. I1=MLENT1.LECT(I1)
  206. I2=MLENT2.LECT(I2)
  207. I3=MLENT1.LECT(I3)
  208. V=(MPOVA1.VPOCHA(I1,1)+MPOVA1.VPOCHA(I3,1) )*0.5D0
  209. VPOCHA(I2,1)=V
  210. 2 CONTINUE
  211.  
  212. SEGDES MELEME
  213. SEGDES MPOVA1,MPOVAL
  214. SEGSUP MLENT1,MLENT2
  215.  
  216. ELSEIF(IKAS.EQ.3)THEN
  217. C CENTRE
  218. CALL LEKTAB(MTABD,'XXVOLUM',MCHPOI)
  219. IF(MCHPOI.EQ.0)RETURN
  220.  
  221. ELSEIF(IKAS.EQ.5.OR.IKAS.EQ.6)THEN
  222. C CENTREP1 et CENTREP0
  223.  
  224. CALL ACME(MTABD,'INEFMD',INEFMD)
  225.  
  226. IF(INEFMD.EQ.1)THEN
  227. C% Le type d'élément fini %m1:8 ne convient pas.
  228. MOTERR( 1: 8) = 'LINE'
  229. CALL ERREUR(927)
  230. RETURN
  231. ENDIF
  232.  
  233. TYPE = ' '
  234. CALL ACMO(MTABD,'XXPSOML',TYPE,MCHELM)
  235. TYPE = ' '
  236. IF (TYPE.NE.'MCHAML ') THEN
  237. CALL ECROBJ('TABLE',MTABD)
  238. CALL KPSOML
  239. TYPE = 'MCHAML'
  240. CALL LIROBJ(TYPE,MCHELM,1,IRET)
  241. IF (IRET.EQ.0)THEN
  242. NOM='XXPSOML'
  243. GOTO 5000
  244. ENDIF
  245. CALL ECMO(MTABD,'XXPSOML','MCHAML',MCHELM)
  246. ENDIF
  247. SEGACT MCHELM
  248.  
  249. TYPE=' '
  250. IF(INEFMD.EQ.2)THEN
  251. CALL ACMO(MTABD,'MACRO1 ',TYPE,MELEME)
  252. IF(TYPE.NE.'MAILLAGE')RETURN
  253. ELSEIF(INEFMD.EQ.3)THEN
  254. CALL ACMO(MTABD,'MAILLAGE',TYPE,MELEME)
  255. NOM='MAILLAGE'
  256. IF(TYPE.NE.'MAILLAGE')GO TO 5000
  257. ENDIF
  258.  
  259. TYPE=' '
  260. CALL ACMO(MTABD,'CENTREP1',TYPE,MELEMS)
  261. IF(TYPE.NE.'MAILLAGE')THEN
  262. CALL KCTRP1(MTABD,MELEMS,1)
  263. ENDIF
  264. CALL KRIPAD(MELEMS,MLENTI)
  265. CALL ACMO(MTABD,'ELTP1NC ',TYPE,MELEMP)
  266. NOM='ELTP1NC '
  267. IF(TYPE.NE.'MAILLAGE')GO TO 5000
  268. CALL CRCHPT('CENTREP1',MELEMS,1,MCHPOI)
  269.  
  270.  
  271. CALL LICHTM(MCHPOI,MPOVAL,TYPC,IGEOM)
  272. SEGACT MELEME
  273. NBSOUS=LISOUS(/1)
  274. IF(NBSOUS.EQ.0)NBSOUS=1
  275. NUTOEL=0
  276.  
  277. NPTD=VPOCHA(/1)
  278. IES=IDIM
  279. MP10=0
  280.  
  281. DO 11 L=1,NBSOUS
  282. IPT1=MELEME
  283. IF(NBSOUS.NE.1)IPT1=LISOUS(L)
  284. SEGACT IPT1
  285.  
  286. MCHAML=ICHAML(L)
  287. SEGACT MCHAML
  288. MELVAL=IELVAL(1)
  289. SEGACT MELVAL
  290.  
  291. NP =IPT1.NUM(/1)
  292. NBEL=IPT1.NUM(/2)
  293.  
  294. IPT2=MELEMP
  295. IF(NBSOUS.NE.1)IPT2=LISOUS(L)
  296. SEGACT IPT2
  297.  
  298. IF(INEFMD.EQ.2.AND.IKAS.EQ.5)NOM0=NOMS(IPT1.ITYPEL)//'MCP1'
  299. IF(INEFMD.EQ.3.AND.IKAS.EQ.5)NOM0=NOMS(IPT1.ITYPEL)//'PRP1'
  300.  
  301. CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  302.  
  303.  
  304. SEGACT IZFFM*MOD
  305. IZHR=KZHR(1)
  306. IZH2=KZHR(2)
  307. SEGACT IZHR*MOD,IZH2*MOD
  308. NES=GR(/1)
  309. NPG=GR(/3)
  310. IZF1=KTP(1)
  311. SEGACT IZF1*MOD
  312. MP1=IZF1.FN(/1)
  313. NP = IPT1.NUM(/1)
  314.  
  315. DO 21 K=1,NBEL
  316. DO 109 I=1,NP
  317. J=IPT1.NUM(I,K)
  318. JC = (J-1)*(IDIM+1)
  319. DO 109 N=1,IDIM
  320. XYZ(N,I)=XCOOR( JC + N )
  321. 109 CONTINUE
  322.  
  323. CALL CALJBC(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE)
  324.  
  325. DO 39 M=1,MP1
  326. M11=LECT(IPT2.NUM(M,K))
  327. M1=M+MP10
  328. c IF(KPOIND.EQ.5)M1=M11
  329.  
  330. U=0.D0
  331. DO 33 LL=1,NPG
  332. U=U+IZF1.FN(M,LL)*PGSQ(LL)
  333. 33 CONTINUE
  334.  
  335. VPOCHA(M1,1)=VPOCHA(M1,1)+U
  336. 39 CONTINUE
  337.  
  338. MP10=MP10+MP1
  339. 21 CONTINUE
  340.  
  341. SEGSUP IZFFM,IZF1,IZHR,IZH2
  342.  
  343.  
  344. SEGDES IPT1,MCHAML,MELVAL,IPT2
  345. NUTOEL=NUTOEL+NBEL
  346.  
  347. 11 CONTINUE
  348.  
  349. SEGSUP MLENTI
  350.  
  351. ENDIF
  352.  
  353.  
  354. SEGDES MTABD,MPOVAL,MCHPOI
  355. CALL ECROBJ('CHPOINT',MCHPOI)
  356. RETURN
  357. 5000 CONTINUE
  358. C Indice %m1:8 : Problème de données détecté dans lektab
  359. IPOINT = 0
  360. MOTERR(1:8) = NOM
  361. CALL ERREUR(792)
  362. RETURN
  363. 1001 FORMAT(20(1X,I5))
  364. 1002 FORMAT(10(1X,1PE11.4))
  365. 1003 FORMAT(6(1X,I7,1X,1PE11.4))
  366. END
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  

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