Télécharger la plaquette Cast3M


fr en

Télécharger gaunew.eso

Retour à la liste

Numérotation des lignes : 

gaunew

  1. C GAUNEW    SOURCE    CB215821  20/11/25    13:29:43     10792          
  2.       SUBROUTINE GAUNEW
  3. C
  4. C     Determination de la direction de descente pour un probleme
  5. C        de moindre carre par la methode de Gauss-Newton
  6. C
  7. C     Possiblilite d'introduire un parametre de viscosite (facultatif)
  8. C
  9. C   TAB2=GANE TAB1 ('AMOR' FLOT1);
  10. C
  11. C   CHPO1 RIGI1=GANE TAB1 'MATR' ('AMOR' FLOT1);
  12. C
  13.       IMPLICIT INTEGER(I-N)
  14.       IMPLICIT REAL*8 (A-H,O-Z)
  15.  
  16. -INC PPARAM
  17. -INC CCOPTIO
  18. -INC SMRIGID
  19. -INC SMCHPOI
  20. -INC SMELEME
  21. -INC CCHAMP
  22. -INC SMLENTI
  23. -INC SMLREEL
  24. -INC SMCOORD
  25. C
  26.       CHARACTER*8 TYPOBJ,CHARR0
  27.       CHARACTER*30 CHARRE
  28.       REAL*8 X0,X1
  29.       LOGICAL L0,L1
  30. *
  31.       PARAMETER(NOPTIO=2)
  32.       CHARACTER*4 OPTION(NOPTIO)
  33.       LOGICAL LMATR
  34.       DATA OPTION/'AMOR','MATR'/
  35. C
  36. C     LECTURE DES OPTIONS
  37. C
  38.       LMATR=.FALSE.
  39.  1    CALL LIRMOT(OPTION,NOPTIO,IRETOU,0)
  40.       IF(IRETOU.EQ.0)GOTO 9
  41.       GOTO(2,3),IRETOU
  42. C
  43. C     OPTION AMOR: On introduit un amortissement
  44. C
  45.  2    CALL LIRREE(AMOR,1,IRETOU)
  46.       IF(IERR.NE.0)RETURN
  47.       GOTO 1
  48. C
  49. C     OPTION MATR: Osort avec la matrice et le second membre
  50. C
  51.  3    LMATR=.TRUE.
  52.       GOTO 1
  53. C
  54.  9    CONTINUE
  55. C
  56. C     LECTURE DE LA TABLE
  57. C
  58.       CALL LIROBJ('TABLE',MTABLE,1,IRETOU)
  59.       IF(IERR.NE.0)RETURN
  60. C
  61. C     INDEX TAB1.SOUSTYPE
  62. C
  63.       CALL ACCTAB(MTABLE,'MOT',I0,X0,'SOUSTYPE',L0,IPO0,
  64.      $                   'MOT',I1,X1,CHARRE,L1,IPO1)
  65.       IF(IERR.NE.0) THEN
  66.         IRET=1
  67.         RETURN
  68.       ENDIF
  69.       IF(CHARRE(1:7).NE.'VECTEUR')THEN
  70.         CALL ERREUR(648)
  71.         IRET=1
  72.         RETURN
  73.       ENDIF
  74. C
  75. C     CHPOINT INDEX 0 (valeur de la fonction)
  76. C
  77.       CALL ACCTAB(MTABLE,'ENTIER',0,X0,'SOUSTYPE',L0,IPO0,
  78.      $                   'CHPOINT',I1,X1,CHARRE,L1,MCHPOI)
  79.       IF(IERR.NE.0) THEN
  80.         IRET=1
  81.         RETURN
  82.       ENDIF
  83. C
  84. C     VERIFICATION DU CARACTERE ELEMENTAIRE DU CHPOINT,
  85. C     SANS SERIE DE FOURIER
  86. C
  87.       SEGACT,MCHPOI
  88.       IF(IPCHP(/1).NE.1)THEN
  89.         WRITE(IOIMP,*)'GAuss NEwton: work only with elementary CHPO'
  90.         MOTERR(1:8)='CHPOINT '
  91.         INTERR(1)=MCHPOI
  92.         CALL ERREUR(110)
  93.         GOTO 9999
  94.       ENDIF
  95.       IF(IFOPOI.EQ.1)THEN
  96.         MOTERR(1:8)='GANE    '
  97.         CALL ERREUR(333)
  98.         GOTO 9999
  99.       ENDIF
  100. C
  101. C     ACTIVATIONS DIVERSES pour la fonction
  102. C
  103.       MSOUPO=IPCHP(1)
  104.       SEGACT,MSOUPO
  105.       MELEME=IGEOC
  106.       MPOVAL=IPOVAL
  107.       SEGACT,MELEME,MPOVAL
  108.       N1=VPOCHA(/1)
  109.       NC1=VPOCHA(/2)
  110. C
  111. C     ON PASSE AUX AUTRES INDEX ENTIERS DE TYPE CHPO CONSECUTIFS DE LA TABLE
  112. C
  113.       JG=0
  114.       SEGINI,MLENTI
  115. 10    CONTINUE
  116.         TYPOBJ=' '
  117.         JG=JG+1
  118.         CALL ACCTAB(MTABLE,'ENTIER',JG,X0,CHARR0,L0,IPO0,
  119.      $                       TYPOBJ,I1,X1,CHARRE,L1,IPO1)
  120.         IF(TYPOBJ.EQ.'CHPOINT')THEN
  121.           SEGADJ,MLENTI
  122.           LECT(JG)=IPO1
  123.           GOTO 10
  124.         ENDIF
  125.       CONTINUE
  126. C
  127. C     ... ET ON VERIFIE QUE TOUT BAIGNE
  128. C
  129.       JG=LECT(/1)
  130.       IF(JG.EQ.0)THEN
  131.         WRITE(IOIMP,*)'GAuss NEwton: The gradient is empty'
  132.         CALL ERREUR(314)
  133.         GOTO 9998
  134.       ENDIF
  135. C
  136.       DO IE1=1,JG
  137.         MCHPO1=LECT(IE1)
  138. C
  139. C     --> VERIFICATION DU CARACTERE ELEMENTAIRE DU CHPOINT,
  140. C         SANS SERIE DE FOURIER
  141. C
  142.         SEGACT,MCHPO1
  143.         IF(MCHPO1.IPCHP(/1).NE.1)THEN
  144.           WRITE(IOIMP,*)'GAuss NEwton: work only with elementary CHPO'
  145.           MOTERR(1:8)='CHPOINT '
  146.           INTERR(1)=MCHPOI
  147.           CALL ERREUR(110)
  148.           SEGDES,MCHPO1
  149.           GOTO 9998
  150.         ENDIF
  151.         IF(MCHPO1.IFOPOI.EQ.1)THEN
  152.           MOTERR(1:8)='GANE    '
  153.           CALL ERREUR(333)
  154.           SEGDES,MCHPO1
  155.           GOTO 9998
  156.         ENDIF
  157. C
  158. C     --> MEME COMPOSANTES, MEME MAILLAGE
  159. C
  160.         MSOUP1=MCHPO1.IPCHP(1)
  161.         SEGACT,MSOUP1
  162. C
  163.         IF (NC1.NE.MSOUP1.NOHARM(/1))THEN
  164.           WRITE(IOIMP,*)
  165.      >    'GAuss NEwton: all CHPO should have the same nb of components'
  166.           CALL ERREUR(665)
  167.           SEGDES,MSOUP1,MCHPO1
  168.           GOTO 9998
  169.         ENDIF
  170. C
  171.         DO IE2=1,NC1
  172.           IF(NOCOMP(IE2).NE.MSOUP1.NOCOMP(IE2))THEN
  173.           WRITE(IOIMP,*)
  174.      >    'GAuss NEwton: all CHPO should have the same components',
  175.      >    ' in the same order'
  176.             CALL ERREUR(488)
  177.             SEGDES,MSOUP1,MCHPO1
  178.             GOTO 9998
  179.           ENDIF
  180.         ENDDO
  181. C
  182.         IPT1=MSOUP1.IGEOC
  183.         SEGACT,IPT1
  184.         IF(N1.NE.IPT1.ICOLOR(/1))THEN
  185.           WRITE(IOIMP,*)
  186.      >    'GAuss NEwton: all CHPO should have the same number',
  187.      >    ' of points'
  188.           MOTERR(1:8) ='MAILLAGE'
  189.           MOTERR(9:16)='NOEUD    '
  190.           CALL ERREUR(403)
  191.           SEGDES,MSOUP1,MCHPO1,IPT1
  192.           GOTO 9998
  193.         ENDIF
  194. C
  195.         DO IE2=1,N1
  196.           IF(NUM(1,IE2).NE.IPT1.NUM(1,IE2))THEN
  197.             WRITE(IOIMP,*)
  198.      >      'GAuss NEwton: all CHPO should have the same',
  199.      >      ' points in the same order'
  200.             CALL ERREUR(403)
  201.             SEGDES,MSOUP1,MCHPO1,IPT1
  202.             GOTO 9998
  203.           ENDIF
  204.         ENDDO
  205. C
  206.         LECT(IE1)=MSOUP1.IPOVAL
  207.         SEGDES,MSOUP1,MCHPO1,IPT1
  208. C
  209.       ENDDO
  210. C
  211. C     ON PREPARE LA MATRICE ET LE SECOND MEMBRE
  212. C
  213. C     1) NOEUDS DU SUPER ELEMENT ET DU CHPO
  214. C
  215.       segact mcoord*mod
  216.       NBPTT=nbpts
  217.       NBPTS=NBPTT+JG
  218.       SEGADJ,MCOORD
  219. C
  220. C     2) SUPER ELEMENT
  221. C
  222.       NBSOUS=0
  223.       NBELEM=1
  224.       NBNN=JG
  225.       NBREF=0
  226.       SEGINI,IPT1
  227.       IPT1.ITYPEL=28
  228.       DO IE1=1,JG
  229.         IPT1.NUM(IE1,1)=NBPTT+ IE1
  230.       ENDDO
  231.       SEGDES,IPT1
  232. C
  233. C     4) DECRIPTEUR POUR LA RIGIDITE
  234. C
  235.       NLIGRP=JG
  236.       NLIGRD=JG
  237.       SEGINI,DESCR
  238.       DO IE1=1,JG
  239.         LISINC(IE1)=NOMDD(11)
  240.         LISDUA(IE1)=NOMDU(11)
  241.         NOELEP(IE1)=IE1
  242.         NOELED(IE1)=IE1
  243.       ENDDO
  244.       SEGDES,DESCR
  245. C
  246. C     5) CONTENU DE LA RIGIDITE
  247. C
  248.       nelrig=1
  249.       SEGINI,XMATRI
  250.       DO IE1=1,NLIGRP
  251.         MPOVA1=LECT(IE1)
  252.         SEGACT,MPOVA1
  253.         DO IE2=1,IE1
  254.           MPOVA2=LECT(IE2)
  255.           SEGACT,MPOVA2
  256. C
  257.           XVAL=0.D0
  258.           DO IE3=1,NC1
  259.             DO IE4=1,N1
  260.               XVAL=XVAL+MPOVA1.VPOCHA(IE4,IE3)*MPOVA2.VPOCHA(IE4,IE3)
  261.             ENDDO
  262.           ENDDO
  263.           RE(IE1,IE2,1)=XVAL
  264.           RE(IE2,IE1,1)=XVAL
  265. C
  266.           SEGDES,MPOVA2
  267.         ENDDO
  268.       ENDDO
  269. C
  270. C     6) CORRECTION EVENTUELLE DE RIGIDITE SUR LA DIAGONALE
  271. C
  272.       IF(AMOR.NE.0.D0)THEN
  273.         DO IE1=1,NLIGRP
  274.           RE(IE1,IE1,1)=RE(IE1,IE1,1)+AMOR
  275.         ENDDO
  276.       ENDIF
  277.       SEGDES,XMATRI
  278. C
  279. C     7) IMATRI
  280. C
  281. *      NELRIG=1
  282. *      SEGINI,IMATRI
  283. *      IMATTT(1)=XMATRI
  284. *      SEGDES,IMATRI
  285. C
  286. C     8) CHAPEAU MRIGID DE LA RIGIDITE
  287. C
  288.       NRIGEL=1
  289.       NRIGE=8
  290.       SEGINI,MRIGID
  291.       MTYMAT='RIGIGITE'
  292.       COERIG(1)=1.D0
  293.       IRIGEL(1,1)=IPT1
  294.       IRIGEL(2,1)=0
  295.       IRIGEL(3,1)=DESCR
  296.       IRIGEL(4,1)=xMATRI
  297.       IRIGEL(5,1)=0
  298.       IRIGEL(6,1)=0
  299.       IRIGEL(7,1)=0
  300.       IRIGEL(8,1)=0
  301.       ICHOLE=0
  302.       IMGEO1=0
  303.       IMGEO2=0
  304.       IFORIG=IFOPOI
  305.       ISUPEQ=0
  306.       SEGDES,MRIGID
  307. C
  308. C     9) SUPPORT DU CHPO
  309. C
  310.       NBSOUS=0
  311.       NBELEM=JG
  312.       NBNN=1
  313.       NBREF=0
  314.       SEGINI,IPT1
  315.       IPT1.ITYPEL=1
  316.       DO IE1=1,JG
  317.         IPT1.NUM(1,IE1)=NBPTT+ IE1
  318.       ENDDO
  319.       SEGDES,IPT1
  320. C
  321. C     10) VALEUR DU CHPO
  322. C
  323.       N=JG
  324.       NC=1
  325.       SEGINI,MPOVA2
  326. C
  327.       DO IE1=1,N
  328.         MPOVA1=LECT(IE1)
  329.         SEGACT,MPOVA1
  330. C
  331.         XVAL=0.D0
  332.         DO IE2=1,NC1
  333.           DO IE3=1,N1
  334.             XVAL=XVAL-VPOCHA(IE3,IE2)*MPOVA1.VPOCHA(IE3,IE2)
  335.           ENDDO
  336.         ENDDO
  337.         MPOVA2.VPOCHA(IE1,1)=XVAL
  338. C
  339.         SEGDES,MPOVA1
  340.       ENDDO
  341. C
  342.       SEGDES,MPOVA2
  343. C
  344. C     11) MSOUPO
  345. C
  346.       SEGINI,MSOUP1
  347.       MSOUP1.NOCOMP(1)=NOMDU(11)
  348.       MSOUP1.IGEOC=IPT1
  349.       MSOUP1.IPOVAL=MPOVA2
  350.       SEGDES,MSOUP1
  351. C
  352. C     12) MCHPOI
  353. C
  354.       NAT=2
  355.       NSOUPO=1
  356.       SEGINI,MCHPO1
  357.       MCHPO1.JATTRI(1)=JATTRI(1)
  358.       MCHPO1.IFOPOI=IFOPOI
  359.       MCHPO1.IPCHP(1)=MSOUP1
  360.       SEGDES,MCHPO1
  361. C
  362. C     RETOUR A GIBIANE
  363. C
  364.       CALL ECROBJ('RIGIDITE',MRIGID)
  365.       CALL ECROBJ('CHPOINT ',MCHPO1)
  366. C
  367. C     ON SORT SI LAMATR
  368. C
  369.       IF(LMATR)GOTO 9998
  370. C
  371. C     DESACTIVATIONS
  372. C
  373.       SEGDES,MELEME,MPOVAL,MSOUPO
  374.       SEGSUP,MLENTI
  375.       SEGDES,MCHPOI
  376. C
  377. C     ON APPELLE RESO
  378. C
  379.       CALL RESOU
  380. C
  381. C     ... ET ON REGARDE LE RESULTAT
  382. C
  383.       CALL LIROBJ('CHPOINT',MCHPOI,1,IRETOU)
  384.       IF(IERR.NE.0)RETURN
  385. C
  386.       SEGACT,MCHPOI
  387.       MSOUPO=IPCHP(1)
  388.       SEGACT,MSOUPO
  389.       MELEME=IGEOC
  390.       MPOVAL=IPOVAL
  391.       SEGACT,MELEME,MPOVAL
  392. C
  393.       SEGINI,MLENTI,MLENT1
  394.       DO IE1=1,JG
  395.         LECT(IE1)=NUM(1,IE1)
  396.         MLENT1.LECT(IE1)=IE1
  397.       ENDDO
  398.       CALL GENOR2(LECT(1),MLENT1.LECT(1),JG)
  399. C
  400.       SEGINI,MLREEL
  401.       DO IE1=1,JG
  402.         PROG(IE1)=VPOCHA(MLENT1.LECT(IE1),1)
  403.       ENDDO
  404.       CALL VECTAB(MLREEL,JG,MTABLE)
  405. C
  406.       CALL ECROBJ('TABLE   ',MTABLE)
  407. C
  408. C     DESTRUCTION
  409. C
  410.       SEGSUP,MLENTI,MLENT1,MLREEL
  411. C
  412.       SEGSUP,MELEME,MPOVAL,MSOUPO,MCHPOI
  413. C
  414.       MCHPOI=MCHPO1
  415.       SEGACT,MCHPOI
  416.       MSOUPO=IPCHP(1)
  417.       SEGACT,MSOUPO
  418.       MELEME=IGEOC
  419.       MPOVAL=IPOVAL
  420.       SEGSUP,MELEME,MPOVAL,MSOUPO,MCHPOI
  421. C
  422.       SEGACT,MRIGID
  423.       MELEME=IRIGEL(1,1)
  424.       DESCR=IRIGEL(3,1)
  425.       xMATRI=IRIGEL(4,1)
  426. *      SEGACT,IMATRI
  427. *      XMATRI=IMATTT(1)
  428.       SEGSUP,MELEME,DESCR,XMATRI,MRIGID
  429. C
  430.       NBPTS=NBPTT
  431.       SEGADJ,MCOORD
  432. C
  433.       RETURN
  434. C
  435. C     ERREURS
  436. C
  437.  9998 SEGDES,MELEME,MPOVAL,MSOUPO
  438.       SEGSUP,MLENTI
  439.  9999 SEGDES,MCHPOI
  440.       RETURN
  441.       END
  442.  
  443.  
  444.  
  445.  
  446.  
  447.  
  448.  
  449.  
  450.  
  451.  
  452.  
  453.  
  454.  
  455.  

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