Télécharger la plaquette Cast3M


fr en

Télécharger knol.eso

Retour à la liste

Numérotation des lignes : 

knol

  1. C KNOL      SOURCE    CB215821  24/04/12    21:16:32     11897          
  2.       SUBROUTINE KNOL
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5. C*************************************************************************
  6. C     Operateur KNOL
  7. C
  8. C     Objet Transforme un CHAMPOINT SOMMET en un CHAMPOINT CENTRE
  9. C
  10. C     SYNTAXE : CHPC = KNOL TABDOM CHPS ;
  11. C
  12. C     TABDOM : Table DOMAINE contenant les supports geometriques de CHPS
  13. C     CHPS   : CHAMPOINT SOMMET
  14. C     CHPC   : CHAMPOINT CENTRE
  15. C
  16. C
  17. C*************************************************************************
  18. -INC PPARAM
  19. -INC CCOPTIO
  20. -INC CCREEL
  21. -INC CCGEOME
  22. -INC SMCOORD
  23. -INC SMELEME
  24.       POINTEUR MELEMS.MELEME,MELEMD.MELEME,SPGD.MELEME
  25. -INC SMMODEL
  26. -INC SMCHPOI
  27.       POINTEUR IZB.MCHPOI,IZBB.MPOVAL
  28. -INC SMLENTI
  29. -INC SIZFFB
  30.       POINTEUR IZF1.IZFFM,IZH2.IZHR,IZW.IZFFM,IZWH.IZHR
  31.       SEGMENT SAJT
  32.       REAL*8 AJT(IDIM,IDIM,NPG),RF1(NP,MP,IDIM),SM1(NP,IDIM)
  33. c     REAL*8 TN1(NP,IDIM),TN2(NP,IDIM)
  34.       ENDSEGMENT
  35.  
  36.       PARAMETER (NTO=4,NBMO=4)
  37.       DIMENSION ITABO(NTO)
  38.       CHARACTER*4 NOMD4
  39.       CHARACTER*8 TYPE,TYPC,LISMO(NBMO),TYPSPG,MTERR,NOM0
  40.       DATA LISMO/'CENTRE  ','        ', 'CENTREP1','MSOMMET'/
  41. C***
  42.  
  43.       TYPSPG='CENTRE  '
  44.       CALL INITI(ITABO,NTO,0)
  45.  4    CONTINUE
  46.  
  47.       CALL QUETYP(TYPE,0,IRET)
  48. c     write(6,*)' Iret de Quetyp= ',iret,' TYPSPG=',typspg
  49.  
  50.       IF(IRET.EQ.0)THEN
  51.        IF(TYPSPG.EQ.'MSOMMET '.OR.TYPSPG.EQ.'CENTRE  '.OR.
  52.      &    TYPSPG.EQ.'CENTREP1'.OR.TYPSPG.EQ.'        ')THEN
  53.         IF(ITABO(1).EQ.1.AND.ITABO(2).EQ.1)GO TO 52
  54.  
  55.          IF(ITABO(1).NE.1)THEN
  56. C%    Il faut spécifier un objet de type %m1:8 et de sous type %m9:16
  57.          MOTERR(1: 8) = 'CHPOINT '
  58.          MOTERR(9:16) = 'DIFFUS  '
  59.          CALL ERREUR(79)
  60.          ENDIF
  61.  
  62.          IF(ITABO(2).NE.1)THEN
  63. C%    Il faut spécifier un objet de type %m1:8 et de sous type %m9:16
  64.          MOTERR(1: 8) = 'MMODEL  '
  65.          MOTERR(9:16) = '        '
  66.          CALL ERREUR(79)
  67.          ENDIF
  68.  
  69.         ENDIF
  70.         RETURN
  71.        ENDIF
  72. *
  73. * Lecture du CHPOIN
  74. *
  75.       IF(TYPE.EQ.'CHPOINT')THEN
  76.        IF(ITABO(1).NE.0)THEN
  77. C% On a déja lu un objet de type %m1:8
  78.         MOTERR(1: 8) = 'CHPOINT '
  79.         CALL ERREUR(966)
  80.         RETURN
  81.        ENDIF
  82.       ITABO(1)=1
  83.  
  84.       CALL LIROBJ('CHPOINT ',IZB,1,IRETOU)
  85.       CALL ACTOBJ('CHPOINT ',IZB,1)
  86.       SEGACT IZB
  87.        IF(IZB.IPCHP(/1).NE.1)THEN
  88. C% Erreur dans le partitionnement
  89.         CALL ERREUR(920)
  90.         RETURN
  91.        ENDIF
  92.       GO TO 4
  93.       ENDIF
  94. *
  95. * Lecture de l'objet modele 'Navier-Stokes'
  96. *
  97. C***
  98.       IF(TYPE.EQ.'MMODEL  ')THEN
  99.        IF(ITABO(2).NE.0)THEN
  100. C% On a déja lu un objet de type %m1:8
  101.         MOTERR(1: 8) = 'MMODEL  '
  102.         MOTERR(9:16) = '        '
  103.         CALL ERREUR(966)
  104.         RETURN
  105.        ENDIF
  106.         ITABO(2)=1
  107.  
  108.         CALL LIROBJ('MMODEL  ',MMODEL,1,IRET)
  109.         CALL ACTOBJ('MMODEL  ',MMODEL,1)
  110.         SEGACT MMODEL
  111.         N1=KMODEL(/1)
  112.         DO 41 L=1,N1
  113.         IMODEL=KMODEL(L)
  114.         SEGACT IMODEL
  115.          IF(FORMOD(1).NE.'NAVIER_STOKES')THEN
  116.           IF(FORMOD(1).NE.'DARCY')THEN
  117. C% On veut un modèle de type %m1:16 .
  118.            MOTERR( 1:16) = 'NAVIER_STOKES   '
  119.            CALL ERREUR(719)
  120.            RETURN
  121.           ENDIF
  122.          ENDIF
  123.  41     CONTINUE
  124.  
  125.         CALL LEKMOD(MMODEL,MTABD,INEFMD)
  126. C  /S INEFMD : Type formulation INEFMD=1 LINE,=2 MACRO,=3 QUADRATIQUE
  127. C                               INEFMD=4 LINB
  128. C KPOIN = 0->SOMMET 1-> FACE     2-> CENTRE  3-> CENTREP0 4-> CENTREP1 5-> MSOMMET
  129.       GO TO 4
  130.       ENDIF
  131. *
  132. * Lecture de l'objet table DOMAINE
  133. *
  134. C***
  135.       IF(TYPE.EQ.'TABLE   ')THEN
  136.        IF(ITABO(4).NE.0)THEN
  137. C% On a déja lu un objet de type %m1:8
  138.         MOTERR(1: 8) = 'TABLE   '
  139.         MOTERR(9:16) = '        '
  140.         CALL ERREUR(966)
  141.         RETURN
  142.        ENDIF
  143.         ITABO(2)=1
  144.  
  145.         CALL LIROBJ('TABLE   ',MTABD,1,IRET)
  146.       GO TO 4
  147.       ENDIF
  148. *
  149. * Lecture d'un mot
  150. *
  151.       IF(TYPE.EQ.'MOT     ')THEN
  152.       CALL LIRMOT(LISMO,NBMO,IP,1)
  153.       TYPSPG=LISMO(IP)
  154.       GO TO 4
  155.       ENDIF
  156.  
  157.  
  158.  
  159.  52   CONTINUE
  160. C
  161. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  162. C------------------ Traitement du cas MSOMMET ---------------------------
  163.  
  164.       IF(TYPSPG.EQ.'MSOMMET')THEN
  165.       CALL LICHTL(IZB,IZBB,TYPC,IGEOMB)
  166.       NC=IZBB.VPOCHA(/2)
  167.       TYPE=' '
  168.       CALL ACMO(MTABD,'SOMMET',TYPE,MELEMS)
  169.       IF(TYPE.NE.'MAILLAGE')GO TO 90
  170.  
  171. C
  172. C On verifie Le support du CHP1 (SOMMET)
  173. C
  174.       CALL KRIPAD(IGEOMB,MLENTI)
  175.       CALL VERPAD(MLENTI,MELEMS,IRET)
  176.  
  177.        IF(IRET.NE.0)THEN
  178. C     Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  179.          MOTERR(1: 8) = 'SOMMET'
  180.          MOTERR(9:16) = 'CHPOINT '
  181.          SEGSUP MLENTI
  182.          RETURN
  183.        ENDIF
  184.       SEGSUP MLENTI
  185.  
  186.       TYPE=' '
  187.       CALL ACMO(MTABD,'MSOMMET',TYPE,MELEME)
  188.       IF(TYPE.NE.'MAILLAGE')GO TO 90
  189.  
  190.       CALL ECROBJ('MAILLAGE',MELEME)
  191.       CALL ECROBJ('CHPOINT ',IZB)
  192.       CALL REDU
  193.  
  194.       RETURN
  195.       ENDIF
  196. C-------------- FIN Traitement du cas MSOMMET ---------------------------
  197. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  198.  
  199. C
  200. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  201. C------------------ Traitement du cas CENTRE ----------------------------
  202.  
  203.       IF(TYPSPG.EQ.'CENTRE')THEN
  204.       CALL LICHTL(IZB,IZBB,TYPC,IGEOMB)
  205.       NC=IZBB.VPOCHA(/2)
  206.       TYPE=' '
  207.       CALL ACMO(MTABD,'SOMMET',TYPE,MELEMS)
  208.       IF(TYPE.NE.'MAILLAGE')GO TO 90
  209.  
  210.  
  211.  
  212.       CALL KRIPAD(IGEOMB,MLENTI)
  213.       CALL VERPAD(MLENTI,MELEMS,IRET)
  214.  
  215.        IF(IRET.NE.0)THEN
  216. C     Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  217.          MOTERR(1: 8) = 'SOMMET'
  218.          MOTERR(9:16) = 'CHPOINT '
  219.          SEGSUP MLENTI
  220.          RETURN
  221.        ENDIF
  222.  
  223.       TYPE=' '
  224.       CALL ACMO(MTABD,'CENTRE',TYPE,MELEMC)
  225.       IF(TYPE.NE.'MAILLAGE')GO TO 90
  226.  
  227.       TYPE=' '
  228.       CALL ACMO(MTABD,'MAILLAGE',TYPE,MELEME)
  229.       IF(TYPE.NE.'MAILLAGE')GO TO 90
  230.  
  231.  
  232.       TYPE='CENTRE'
  233. c     write(6,*)' Apparemment on traite le cas centre !!!'
  234.       CALL CRCHPT(TYPE,MELEMC,NC,MCHPOI)
  235.       CALL LICHTM(MCHPOI,MPOVAL,TYPC,IGEOM)
  236.  
  237.       SEGACT MELEME
  238.  
  239.       NBSOUS=LISOUS(/1)
  240.       IF(NBSOUS.EQ.0)NBSOUS=1
  241.  
  242.       KK=0
  243.       DO 1 L=1,NBSOUS
  244.       IPT1=MELEME
  245.       IF(NBSOUS.NE.1)IPT1=LISOUS(L)
  246.       SEGACT IPT1
  247.       NP=IPT1.NUM(/1)
  248.       NEL=IPT1.NUM(/2)
  249.       DO 10 K=1,NEL
  250.       KK=KK+1
  251.       DO 13 N=1,NC
  252.       UU=0.D0
  253.       DO 12 I=1,NP
  254.       IU=LECT(IPT1.NUM(I,K))
  255.       UU=UU+IZBB.VPOCHA(IU,N)
  256.  12   CONTINUE
  257.       VPOCHA(KK,N)=UU/FLOAT(NP)
  258.  13   CONTINUE
  259.  
  260.  10   CONTINUE
  261.  1    CONTINUE
  262. C
  263.       SEGSUP MLENTI
  264.       CALL ECROBJ('CHPOINT ',MCHPOI)
  265.  
  266.       RETURN
  267.       ENDIF
  268. C-------------- FIN Traitement du cas CENTRE ----------------------------
  269. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  270.  
  271. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  272. C------------------ Traitement du cas CENTREP1 --------------------------
  273.  
  274.       IF(TYPSPG.EQ.'CENTREP1')THEN
  275.  
  276.           MTERR='EF CTRP1'
  277.           IF(INEFMD.EQ.2)NOMD4='MCP1'
  278.           IF(INEFMD.EQ.3)NOMD4='PRP1'
  279.           IF(INEFMD.NE.2.AND.INEFMD.NE.3)THEN
  280. C        Option %m1:8 incompatible avec les données
  281.             MOTERR( 1: 8) = MTERR
  282.             CALL ERREUR(803)
  283.             RETURN
  284.           ENDIF
  285.  
  286.           CALL LEKTAB(MTABD,'CENTREP1',SPGD)
  287.           CALL LEKTAB(MTABD,'ELTP1NC ',MELEMD)
  288.  
  289.  
  290.       CALL LICHTL(IZB,IZBB,TYPC,IGEOMB)
  291.       NC=IZBB.VPOCHA(/2)
  292.       TYPE=' '
  293.       CALL ACMO(MTABD,'SOMMET',TYPE,MELEMS)
  294.       IF(TYPE.NE.'MAILLAGE')GO TO 90
  295.  
  296.  
  297.       CALL KRIPAD(IGEOMB,MLENT1)
  298.       CALL VERPAD(MLENT1,MELEMS,IRET)
  299.  
  300.        IF(IRET.NE.0)THEN
  301. C     Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  302.          MOTERR(1: 8) = 'SOMMET'
  303.          MOTERR(9:16) = 'CHPOINT '
  304.          SEGSUP MLENTI
  305.          RETURN
  306.        ENDIF
  307.  
  308.       TYPE=' '
  309.       CALL ACMO(MTABD,'MAILLAGE',TYPE,MELEME)
  310.       IF(TYPE.NE.'MAILLAGE')GO TO 90
  311.  
  312.       IF(INEFMD.EQ.2)CALL LEKTAB(MTABD,'MACRO1',MELEME)
  313.       IF (IERR.NE.0) RETURN
  314.  
  315.       TYPE='CENTREP1'
  316. c     write(6,*)' Apparemment on traite le cas centrep1 !!!'
  317.       CALL CRCHPT(TYPE,SPGD,NC,MCHPOI)
  318.       CALL LICHTM(MCHPOI,MPOVAL,TYPC,IGEOM)
  319. c..........................................................................
  320.  
  321.  
  322. c     IK3=0
  323.       IAXI=0
  324.       IF(IFOMOD.EQ.0)IAXI=2
  325.       DEUPI=1.D0
  326.       IF(IAXI.NE.0)DEUPI=2.D0*XPI
  327.  
  328.       NC=MPOVAL.VPOCHA(/2)
  329.       CALL KRIPAD(SPGD,MLENTI)
  330.  
  331.       SEGACT MELEME
  332.  
  333.       NKD=0
  334.       DO 101 L=1,MAX(1,LISOUS(/1))
  335.       SEGACT MELEMD
  336.       IPT1=MELEME
  337.       IPT2=MELEMD
  338.       IF(LISOUS(/1).NE.0)IPT1=LISOUS(L)
  339.       SEGACT IPT1
  340.       IF(MELEMD.LISOUS(/1).NE.0)THEN
  341.       IPT2=MELEMD.LISOUS(L)
  342.       NKD=0
  343.       ENDIF
  344.       SEGACT IPT2
  345.       MP=IPT2.NUM(/1)
  346.  
  347.       NOM0 = NOMS(IPT1.ITYPEL)//NOMD4
  348. c     write(6,*)' KNOL 1er KALPBG NOM0=',NOM0,IPT1
  349.       CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  350.       IF(IZFFM.EQ.0)RETURN
  351.       SEGACT IZFFM*MOD
  352.       IZHR=KZHR(1)
  353.       SEGACT IZHR*MOD
  354.       IZF1 = KTP(1)
  355.       IZH2 = KZHR(2)
  356.       IZW  = IZF1
  357.       SEGACT IZW*MOD
  358.       IF(MP.NE.IZW.FN(/1))THEN
  359.       write(6,*)' Gross problem ds KNOL '
  360.       write(6,*)' NOM0=',NOM0 ,' NOMD4=',NOMD4
  361.       write(6,*)' MP=',MP,' IZW.FN(/1)='
  362.      & ,IZW.FN(/1)
  363.       return
  364.       ENDIF
  365.  
  366.       NES=GR(/1)
  367.       NPG=GR(/3)
  368.       NP = IPT1.NUM(/1)
  369.       NBEL=IPT1.NUM(/2)
  370.       SEGINI SAJT
  371. c      write(6,*)' AVANT 108 NC=',NC,' NBEL=',NBEL,MP,NP,NC
  372. c      write(6,*)' AVANT 108 IK4=',IK4
  373.         DO 108 KE=1,NBEL
  374.  
  375.         NKD=NKD+1
  376.  
  377.         DO 109 I=1,NP
  378.         J=IPT1.NUM(I,KE)
  379.         DO 109 N=1,IDIM
  380.         XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
  381.  109    CONTINUE
  382.  
  383.         CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,
  384.      *  IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  385.  
  386.         CALL INITD(SM1,(MP*IDIM),0.D0)
  387.  
  388. C=======================================================================
  389.  
  390. C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  391. C...... Source
  392.         DO 710 I=1,MP
  393.           U2=0.D0
  394.           U4=0.D0
  395.          DO 717 N=1,NC
  396.           DO 715 LG=1,NPG
  397.           WT=IZW.FN(I,LG)
  398.           U4=U4+WT*PGSQ(LG)*DEUPI*RPG(LG)
  399.           UJ=0.D0
  400.           DO 714 J=1,NP
  401.           JU=MLENT1.LECT(IPT1.NUM(J,KE))
  402.           UJ=UJ+FN(J,LG)*IZBB.VPOCHA(JU,N)
  403.  714      CONTINUE
  404.           U2=U2+UJ*WT*PGSQ(LG)*DEUPI*RPG(LG)
  405.  715      CONTINUE
  406.  
  407.           SM1(I,N)=SM1(I,N)+(U2/U4)
  408.  717     CONTINUE
  409.  710    CONTINUE
  410. c        write(6,*)' SM1 *****'
  411. c        do 711 n=1,nc
  412. c        write(6,1002)(sm1(i,n),i=1,mp)
  413. c711     continue
  414. C...... Source Fin
  415. C=======================================================================
  416.  
  417.  
  418. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  419. C ...... Chargement Second membre
  420.         DO 910 I=1,MP
  421.          I1=LECT(IPT2.NUM(I,NKD))
  422.         DO 910 N=1,NC
  423.          VPOCHA(I1,N)=VPOCHA(I1,N)+SM1(I,N)
  424.  910    CONTINUE
  425. C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  426.  
  427.  108    CONTINUE
  428.  
  429.  
  430.       SEGSUP IZFFM,IZHR,IZF1,IZH2
  431.       SEGSUP SAJT
  432.  
  433.  101  CONTINUE
  434.       SEGSUP MLENTI
  435.  
  436. c..........................................................................
  437.       SEGSUP MLENTI,MLENT1
  438.       CALL ECROBJ('CHPOINT ',MCHPOI)
  439.  
  440.       RETURN
  441.       ENDIF
  442. C-------------- FIN Traitement du cas CENTREP1 --------------------------
  443. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  444.  
  445.  
  446.  90   CONTINUE
  447.       WRITE(6,*)'Interruption anormale de KLNO '
  448.       RETURN
  449.  
  450.  1001 FORMAT(20(1X,I5))
  451.  1002 FORMAT(10(1X,1PE11.4))
  452.  1008 FORMAT(10(1X,A8))
  453.       END
  454.  
  455.  
  456.  
  457.  
  458.  

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