Télécharger la plaquette Cast3M


fr en

Télécharger chaga1.eso

Retour à la liste

Numérotation des lignes : 

chaga1

  1. C CHAGA1    SOURCE    SP204843  24/03/29    21:15:02     11878          
  2.  
  3. C=======================================================================
  4. C=                            C H A G A 1                              =
  5. C=                            -----------                              =
  6. C=                                                                     =
  7. C=  Fonction :                                                         =
  8. C=  ----------                                                         =
  9. C=                                                                     =
  10. C=  Creation du MCHAML de SOURCE de chaleur GAUSSIENNE decrit par les  =
  11. C=  caracteristiques du modele passees en entree.                      =
  12. C=                                                                     =
  13. C=                                                                     =
  14. C=  Entrees :                                                          =
  15. C=  ---------                                                          =
  16. C=                                                                     =
  17. C=  IPMCHA : pointeur sur segment MCHAML de la sous-zone traitee       =
  18. C=                                                                     =
  19. C=  IPMAIL : pointeur sur segment MELEME du maillage ou creer la       =
  20. C=           source de chaleur                                         =
  21. C=                                                                     =
  22. C=  IPINTE : pointeur sur segment MINTE de l'element finis             =
  23. C=                                                                     =
  24. C=  IK     : indicateur du "type" de distribution gaussienne :         =
  25. C=           IK = 1 : source gaussienne SPHERIQUE                      =
  26. C=           IK = 2 : source gaussienne ELLIPTIQUE                     =
  27. C=           IK = 3 : source gaussienne ELARGIE                        =
  28. C=                                                                     =
  29. C=                                                                     =
  30. C=  Sortie :                                                           =
  31. C=  --------                                                           =
  32. C=                                                                     =
  33. C=  IPQGAU : pointeur sur segment MELVAL des valeurs de la source de   =
  34. C=           chaleur GAUSSIENNE                                        =
  35. C=                                                                     =
  36. C=  XQ0    : valeur de la quantite de chaleur imposee, utilisee en     =
  37. C=           sortie dans chamas pour "renormaliser" le champ a QTOT    =
  38. C=                                                                     =
  39. C=======================================================================
  40.  
  41.       SUBROUTINE CHAGA1(IPMCHA,IPMAIL,IPINTE,IK,IPQGAU,XQ0)
  42.  
  43.       IMPLICIT INTEGER(I-N)
  44.       IMPLICIT REAL*8(A-H,O-Z)
  45.  
  46. -INC PPARAM
  47. -INC CCOPTIO
  48. -INC CCREEL
  49. -INC SMCOORD
  50. -INC SMELEME
  51. -INC SMCHAML 
  52. -INC SMINTE  
  53.  
  54.       REAL*8      XCO1(3),XCD1(3),XCX1(3),XCL1(3),XCU1(3)
  55.       CHARACTER*8 TYPC1
  56.       CHARACTER*4 MOCAR(7)
  57.       CHARACTER*(LOCOMP) NOMC1
  58.  
  59.       DATA MOCAR /'QTOT','ORIG','RGAU','DIRE','ZGAU','DIRL','LGAU'/
  60.  
  61.  
  62.       IPQGAU = 0
  63.  
  64. C==== Recuperation des caracteristiques de la source de chaleur
  65. C
  66. C     Initialisation des valeurs des parametres de la source
  67. C     XQ0 : valeur de QTOT
  68. C     NPO1, NPD1 : numeros points origine & direction
  69. C     XR0,  XZ0  : valeurs de RGAU et ZGAU
  70.       XQ0  = 0.D0
  71.       NPO1 = 0
  72.       NPD1 = 0
  73.       XR0  = 0.D0
  74.       XZ0  = 0.D0
  75.       NPL1 = 0
  76.       XL0  = 0.D0
  77. C
  78. C     NCAR1 : nombre de caracteristiques a lire
  79.       NCAR1  = 3
  80.       IF (IK.EQ.2) NCAR1 = 5
  81.       IF (IK.EQ.3) NCAR1 = 7
  82. C
  83.       MCHAM1=IPMCHA
  84. C     SEGACT,MCHAM1
  85.       NCO1 = MCHAM1.IELVAL(/1)
  86.       DO 100 I=1,NCAR1
  87.         NOMC1 = MOCAR(I)
  88.         CALL PLACE(MCHAM1.NOMCHE,NCO1,IPLAC,NOMC1)
  89. C     write(6,*) 'NONC1 =',NOMC1
  90. C     write(6,*) 'MCHAM1.NOMCHE =',(MCHAM1.NOMCHE(ii),ii=1,NCO1)
  91.         IF (IPLAC.EQ.0) THEN
  92.           GOTO 997
  93.         ELSE
  94.           TYPC1 = MCHAM1.TYPCHE(IPLAC)(1:8)
  95.           MELVA1 = MCHAM1.IELVAL(IPLAC)
  96.           IF (MELVA1.VELCHE(/1).GT.1
  97.      &    .OR.MELVA1.VELCHE(/2).GT.1) GOTO 999
  98. C
  99. C         Lecture QTOT
  100.           IF (I.EQ.1) THEN
  101.             IF (TYPC1.NE.'REAL*8  ') GOTO 998
  102.             SEGACT,MELVA1
  103.             XQ0 =  MELVA1.VELCHE(1,1)
  104. C           SEGDES,MELVA1
  105. C           IF (XQ0.LT.XZPREC) THEN
  106. C             MOTERR(1:40)=NOMC1//'                                    '
  107. C             CALL ERREUR(1048)
  108. C             RETURN
  109. C           ENDIF
  110. C
  111. C         Lecture ORIG
  112.           ELSEIF (I.EQ.2) THEN
  113.             IF (TYPC1.NE.'POINTEUR') GOTO 998
  114.             MELVA1 = MCHAM1.IELVAL(IPLAC)
  115.             SEGACT,MELVA1
  116.             NPO1 =  MELVA1.IELCHE(1,1)
  117. C           SEGDES,MELVA1
  118. C
  119. C         Lecture RGAU
  120.           ELSEIF (I.EQ.3) THEN
  121.             IF (TYPC1.NE.'REAL*8  ') GOTO 998
  122.             MELVA1 = MCHAM1.IELVAL(IPLAC)
  123.             SEGACT,MELVA1
  124.             XR0 = MELVA1.VELCHE(1,1)
  125. C           SEGDES,MELVA1
  126.             IF (ABS(XR0).LT.XZPREC) THEN
  127.               MOTERR(1:40)=NOMC1//'                                    '
  128.               CALL ERREUR(1048)
  129.               RETURN
  130.             ENDIF
  131. C
  132. C         Lecture DIRE
  133.           ELSEIF (I.EQ.4) THEN
  134.             IF (TYPC1.NE.'POINTEUR') GOTO 998
  135.             MELVA1 = MCHAM1.IELVAL(IPLAC)
  136.             SEGACT,MELVA1
  137.             NPD1 =  MELVA1.IELCHE(1,1)
  138. C           SEGDES,MELVA1
  139. C
  140. C         Lecture ZGAU
  141.           ELSEIF (I.EQ.5) THEN
  142.             IF (TYPC1.NE.'REAL*8  ') GOTO 998
  143.             MELVA1 = MCHAM1.IELVAL(IPLAC)
  144.             SEGACT,MELVA1
  145.             XZ0 = MELVA1.VELCHE(1,1)
  146. C           SEGDES,MELVA1
  147.             IF (ABS(XZ0).LT.XZPREC) THEN
  148.               MOTERR(1:40)=NOMC1//'                                    '
  149.               CALL ERREUR(1048)
  150.               RETURN
  151.             ENDIF
  152. C
  153. C         Lecture DIRL
  154.           ELSEIF (I.EQ.6) THEN
  155.             IF (TYPC1.NE.'POINTEUR') GOTO 998
  156.             MELVA1 = MCHAM1.IELVAL(IPLAC)
  157.             SEGACT,MELVA1
  158.             NPL1 =  MELVA1.IELCHE(1,1)
  159. C           SEGDES,MELVA1
  160. C
  161. C         Lecture LGAU
  162.           ELSEIF (I.EQ.7) THEN
  163.             IF (TYPC1.NE.'REAL*8  ') GOTO 998
  164.             MELVA1 = MCHAM1.IELVAL(IPLAC)
  165.             SEGACT,MELVA1
  166.             XL0 = MELVA1.VELCHE(1,1)
  167. C           SEGDES,MELVA1
  168. C
  169.           ENDIF
  170.         ENDIF
  171.  100  CONTINUE
  172. C     SEGDES,MCHAM1
  173. C
  174. C     Remplissage du vecteur de coordonnes de ORIG
  175.       IDIMP1 = IDIM+1
  176.       IF (NPO1.EQ.0.OR.NPO1.GT.nbpts) THEN
  177.         NOMC1='ORIG'
  178.         GOTO 998
  179.       ELSE
  180.         DO 120 ID1=1,IDIM
  181.           XCO1(ID1)=XCOOR((NPO1-1)*IDIMP1+ID1)
  182.  120    CONTINUE
  183.       ENDIF
  184. C
  185. C     Remplissage des vecteurs de coordonnes de DIRE et DIRL
  186. *     On orthogonalise DIRL a DIRE
  187. C     XN1 : norme de DIRE au carre
  188. C     XL1 : norme de DIRL au carre
  189.       XN1 = 1.D0
  190.       XL1 = 1.D0
  191.       IF (IK.EQ.1) THEN
  192.         XZ0 = XR0
  193.         XCD1(1)=1.D0
  194.         IF (IDIM.GT.1) THEN
  195.           DO 130 ID1=2,IDIM
  196.             XCD1(ID1)=0.D0
  197.  130      CONTINUE
  198.         ENDIF
  199.       ELSE IF (IK.EQ.2) THEN
  200.         XN1  = 0.D0
  201.         DO 140 ID1=1,IDIM
  202.           XCD1(ID1)=XCOOR((NPD1-1)*IDIMP1+ID1)
  203.           XN1 = XCD1(ID1) ** 2 + XN1
  204.  140    CONTINUE
  205.         IF (XN1.LT.XZPREC) THEN
  206.           CALL ERREUR(239) 
  207.           RETURN
  208.         ENDIF
  209.       ELSE IF (IK.EQ.3) THEN
  210.         XN1  = 0.D0
  211.         XX1  = 0.D0
  212.         DO 150 ID1=1,IDIM
  213.           XCD1(ID1)=XCOOR((NPD1-1)*IDIMP1+ID1)
  214.           XN1 = XCD1(ID1) ** 2 + XN1
  215.           XCL1(ID1)=XCOOR((NPL1-1)*IDIMP1+ID1)
  216.           XX1 = XX1 + XCD1(ID1) * XCL1(ID1)
  217.  150    CONTINUE
  218.         IF (XN1.LT.XZPREC) THEN
  219.           CALL ERREUR(239) 
  220.           RETURN
  221.         ENDIF
  222.         XL1  = 0.D0
  223.         DO 151 IL1= 1,IDIM
  224.           XCL1(IL1) = XCL1(IL1) - XX1 * XCD1(IL1) / XN1
  225.           XL1 = XCL1(IL1) ** 2 + XL1
  226.  151    CONTINUE
  227.         IF (XL1.LT.(XZPREC*XN1)) THEN
  228.           IF (ABS(XX1).LT.XZPREC) THEN
  229.             CALL ERREUR(239) 
  230.           ELSE
  231.             CALL ERREUR(21) 
  232.           ENDIF
  233.           RETURN
  234.         ENDIF
  235.         XCU1(3)=XCL1(1)*XCD1(2)-XCL1(2)*XCD1(1)/SQRT(XN1)/SQRT(XL1)
  236.         XCU1(1)=XCL1(2)*XCD1(3)-XCL1(3)*XCD1(2)/SQRT(XN1)/SQRT(XL1)
  237.         XCU1(2)=XCL1(3)*XCD1(1)-XCL1(1)*XCD1(3)/SQRT(XN1)/SQRT(XL1)
  238.         XU1 = XCU1(1)**2 + XCU1(2)**2 + XCU1(3)**2
  239.       ELSE
  240.         WRITE(IOIMP,*) ' *** Dans CHAGA1, cas (IK) non prevu'
  241.         CALL ERREUR(5)
  242.         RETURN
  243.       ENDIF
  244. C     write(6,*) 'IK =',IK
  245. C     write(6,*) 'XQ0 =',XQ0
  246. C     write(6,*) 'XR0 =',XR0
  247. C     write(6,*) 'XZ0 =',XZ0
  248. c     write(6,*) 'XL0 =',XL0
  249. C     write(6,*) 'NPO1 =',NPO1
  250. C     write(6,*) 'NPD1 =',NPD1
  251. C     write(6,*) 'XCO1 =',(XCO1(i),i=1,IDIM)
  252. C     write(6,*) 'XN1,XCD1 =',XN1,(XCD1(i),i=1,IDIM)
  253. C     write(6,*) 'XL1,XCL1 =',XL1,(XCL1(i),i=1,IDIM)
  254. C     write(6,*) 'XU1,XCU1 =',XU1,(XCU1(i),i=1,IDIM)
  255.  
  256.  
  257. C==== Construction du MELVAL de la distribution Gaussienne de chaleur
  258.  
  259. C     Activation du maillage
  260.       IPT1 = IPMAIL
  261.       SEGACT,IPT1
  262.       NBPT1 = IPT1.NUM(/1)
  263.       NBEL1 = IPT1.NUM(/2)
  264.  
  265. C     Petite verif. sous-zone elementaire
  266.       NBS1 = IPT1.LISOUS(/1)
  267.       IF (NBS1.NE.0) THEN
  268.         WRITE(6,*) 'Dans CHAGA1 : le maillage a plusieurs sous-zones ?'
  269.         CALL ERREUR(21)
  270.         RETURN
  271.       ENDIF
  272. C
  273. C     Activation du segment MINTE
  274. C     write(6,*) 'IPINTE=',IPINTE
  275.       MINTE1 = IPINTE
  276.       SEGACT,MINTE1
  277.  
  278. C     Creation du MELVAL
  279.       N1PTEL = MINTE1.POIGAU(/1)
  280.       N1EL   = NBEL1
  281.       N2PTEL = 0
  282.       N2EL   = 0
  283.       SEGINI,MELVA2
  284.  
  285. C     Facteur de normalisation de la puissance :
  286.       IF (IDIM.EQ.2.AND.IFOUR.NE.0) THEN
  287.         XQ1    = 4.D0 * XQ0 / (XR0 * XZ0 * XPI)
  288.       ELSEIF (IDIM.EQ.2.AND.IFOUR.EQ.0) THEN
  289. C       write (6,*) ' *** Cas axis' 
  290.         RC2    = SQRT (2.D0)
  291.         RCPI   = SQRT (XPI)
  292.         XIK1   = 0.25D0 * XR0 * XR0
  293.         XIK1   = XIK1 * EXP(-2.D0 * XCO1(1) * XCO1(1) / XR0 / XR0)
  294.         XIK2   = 0.5D0 * XR0 * XCO1(1) * RCPI / RC2
  295.         XIK3   = XR0 * XCO1(1) * RCPI / (RC2 * RC2 * RC2)
  296.         XIK3   = XIK3 * ERF(XCO1(1) * RC2 / XR0)
  297.         XQ1    = XPI * XZ0 * RCPI / RC2 
  298.         XQ1    = XQ1 * (XIK1 + XIK2 + XIK3)
  299.         XQ1    = XQ0 / XQ1
  300.       ELSEIF (IDIM.EQ.3) THEN
  301.         IF (IK.EQ.3) THEN
  302.           XR1  = (SQRT ((XPI**3) / 32.D0)) * XR0 * XR0 * XZ0
  303.           XR1  = XR1 + 0.5D0*XPI*XR0*XZ0*ABS(XL0)
  304.           XQ1  = XQ0 / XR1
  305.         ELSE
  306.           XR1  = XR0 * XR0 * XZ0
  307.           XQ1  = SQRT (32.D0 / (XPI ** 3))
  308.           XQ1  = XQ1 * XQ0 / XR1 
  309.         ENDIF
  310.       ELSE
  311.         GOTO 996
  312.       ENDIF
  313. c     write(6,*) 'XQ1=',XQ1
  314.  
  315. C     Boucle sur les elements du maillage
  316.       DO 200 IEL1=1,NBEL1
  317. C
  318. C       Boucle sur les pts de Gauss
  319. C       et calcul de la fonction Gaussienne
  320. C       RR1 : distance a l'origine ORIG au carre
  321. C       XZ1 : hauteur selon la direction DIRE au carre
  322. C       XR1 : rayon dans le plan orthog. a DIRE au carre
  323.         DO 201 IG1=1,N1PTEL
  324.           RR1 = 0.D0
  325.           XR1 = 0.D0
  326.           XZ1 = 0.D0
  327.           VX1 = 0.D0
  328.           UX1 = 0.D0
  329.           DO 203 ID1=1,IDIM
  330. C           XPG1 : coordonnees ID1 du Pt de Gauss
  331.             XPG1 = 0.D0
  332.             DO 202 JPT1=1,NBPT1
  333.               NPTJM1 = IPT1.NUM(JPT1,IEL1)-1
  334.               XPTJ1  = XCOOR(NPTJM1*IDIMP1+ID1)
  335.               VNJ1   = MINTE1.SHPTOT(1,JPT1,IG1)
  336.               XPG1   = XPG1 + XPTJ1*VNJ1
  337.  202        CONTINUE
  338.             XCX1(ID1) = XPG1 - XCO1(ID1)
  339.             XZ1 = XCX1(ID1) * XCD1(ID1) + XZ1
  340.             IF (IK.EQ.3) THEN
  341.               VX1 = XCX1(ID1) * XCL1(ID1) + VX1
  342.               UX1 = XCX1(ID1) * XCU1(ID1) + UX1
  343.             ELSE
  344.               RR1 = XCX1(ID1) ** 2 + RR1
  345.             ENDIF
  346.  203      CONTINUE
  347.           IF (IK.EQ.1.OR.XZ1.LE.0.) THEN
  348.             XZ1 = XZ1*XZ1/ XN1 
  349.             IF (IK.EQ.3) THEN
  350.               UX1 = UX1*UX1/XU1
  351.               XS1 = UX1/XR0/XR0 + XZ1/XZ0/XZ0
  352.               VX1 = SQRT(VX1*VX1/XL1) - XL0
  353.               IF (VX1.GT.0.D0) THEN
  354.                 XS1 = VX1*VX1/XR0/XR0 + XS1
  355.               ENDIF
  356.             ELSE
  357.               XR1 = RR1 - XZ1
  358.               XS1 = XR1/XR0/XR0
  359.               XS1 = XZ1/XZ0/XZ0 + XS1
  360.             ENDIF
  361.             XS1 = EXP(-2.D0 * XS1)
  362.             MELVA2.VELCHE(IG1,IEL1) = XS1 * XQ1
  363.           ELSE
  364.             MELVA2.VELCHE(IG1,IEL1) = 0.D0
  365.           ENDIF
  366.  201    CONTINUE
  367.  200  CONTINUE
  368. C     SEGDES,IPT1,MINTE1,MELVA2
  369. C
  370. C==== Affectation du resultat et sortie
  371.       IPQGAU = MELVA2
  372. C
  373.       RETURN
  374.  
  375. C==== Gestion erreurs et fin
  376.  
  377. C     Option indisponible
  378.  996  CONTINUE
  379.       INTERR(1) = IDIM
  380.       CALL ERREUR(709)
  381.       RETURN
  382.  
  383. C     Il manque la donnee d'une composante
  384.  997  CONTINUE
  385.       CALL ERREUR(80)
  386.       RETURN
  387.  
  388. C     Le type de la composante NOMC1 n'est pas celui attendu
  389.  998  CONTINUE
  390.       MOTERR(1:8)=NOMC1
  391.       CALL ERREUR(679)
  392.       RETURN
  393.  
  394. C     La composante est attendue constante par sous-zones
  395.  999  CONTINUE
  396.       MOTERR(1:4)=NOMC1
  397.       MOTERR(5:20)='CARACTERISTIQUES'
  398.       CALL ERREUR(1065)
  399.       RETURN
  400.  
  401.       END
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  

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