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    PV090527  26/04/28    21:15:36     12529          
  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.       RIGREL=0
  250.       SEGINI,XMATRI
  251.       DO IE1=1,NLIGRP
  252.         MPOVA1=LECT(IE1)
  253.         SEGACT,MPOVA1
  254.         DO IE2=1,IE1
  255.           MPOVA2=LECT(IE2)
  256.           SEGACT,MPOVA2
  257. C
  258.           XVAL=0.D0
  259.           DO IE3=1,NC1
  260.             DO IE4=1,N1
  261.               XVAL=XVAL+MPOVA1.VPOCHA(IE4,IE3)*MPOVA2.VPOCHA(IE4,IE3)
  262.             ENDDO
  263.           ENDDO
  264.           RE(IE1,IE2,1)=XVAL
  265.           RE(IE2,IE1,1)=XVAL
  266. C
  267.           SEGDES,MPOVA2
  268.         ENDDO
  269.       ENDDO
  270. C
  271. C     6) CORRECTION EVENTUELLE DE RIGIDITE SUR LA DIAGONALE
  272. C
  273.       IF(AMOR.NE.0.D0)THEN
  274.         DO IE1=1,NLIGRP
  275.           RE(IE1,IE1,1)=RE(IE1,IE1,1)+AMOR
  276.         ENDDO
  277.       ENDIF
  278.       SEGDES,XMATRI
  279. C
  280. C     7) IMATRI
  281. C
  282. *      NELRIG=1
  283. *      SEGINI,IMATRI
  284. *      IMATTT(1)=XMATRI
  285. *      SEGDES,IMATRI
  286. C
  287. C     8) CHAPEAU MRIGID DE LA RIGIDITE
  288. C
  289.       NRIGEL=1
  290.       NRIGE=8
  291.       SEGINI,MRIGID
  292.       MTYMAT='RIGIGITE'
  293.       COERIG(1)=1.D0
  294.       IRIGEL(1,1)=IPT1
  295.       IRIGEL(2,1)=0
  296.       IRIGEL(3,1)=DESCR
  297.       IRIGEL(4,1)=xMATRI
  298.       IRIGEL(5,1)=0
  299.       IRIGEL(6,1)=0
  300.       IRIGEL(7,1)=0
  301.       IRIGEL(8,1)=0
  302.       ICHOLE=0
  303.       IMGEO1=0
  304.       IMGEO2=0
  305.       IFORIG=IFOPOI
  306.       ISUPEQ=0
  307.       SEGDES,MRIGID
  308. C
  309. C     9) SUPPORT DU CHPO
  310. C
  311.       NBSOUS=0
  312.       NBELEM=JG
  313.       NBNN=1
  314.       NBREF=0
  315.       SEGINI,IPT1
  316.       IPT1.ITYPEL=1
  317.       DO IE1=1,JG
  318.         IPT1.NUM(1,IE1)=NBPTT+ IE1
  319.       ENDDO
  320.       SEGDES,IPT1
  321. C
  322. C     10) VALEUR DU CHPO
  323. C
  324.       N=JG
  325.       NC=1
  326.       SEGINI,MPOVA2
  327. C
  328.       DO IE1=1,N
  329.         MPOVA1=LECT(IE1)
  330.         SEGACT,MPOVA1
  331. C
  332.         XVAL=0.D0
  333.         DO IE2=1,NC1
  334.           DO IE3=1,N1
  335.             XVAL=XVAL-VPOCHA(IE3,IE2)*MPOVA1.VPOCHA(IE3,IE2)
  336.           ENDDO
  337.         ENDDO
  338.         MPOVA2.VPOCHA(IE1,1)=XVAL
  339. C
  340.         SEGDES,MPOVA1
  341.       ENDDO
  342. C
  343.       SEGDES,MPOVA2
  344. C
  345. C     11) MSOUPO
  346. C
  347.       SEGINI,MSOUP1
  348.       MSOUP1.NOCOMP(1)=NOMDU(11)
  349.       MSOUP1.IGEOC=IPT1
  350.       MSOUP1.IPOVAL=MPOVA2
  351.       SEGDES,MSOUP1
  352. C
  353. C     12) MCHPOI
  354. C
  355.       NAT=2
  356.       NSOUPO=1
  357.       SEGINI,MCHPO1
  358.       MCHPO1.JATTRI(1)=JATTRI(1)
  359.       MCHPO1.IFOPOI=IFOPOI
  360.       MCHPO1.IPCHP(1)=MSOUP1
  361.       SEGDES,MCHPO1
  362. C
  363. C     RETOUR A GIBIANE
  364. C
  365.       CALL ECROBJ('RIGIDITE',MRIGID)
  366.       CALL ECROBJ('CHPOINT ',MCHPO1)
  367. C
  368. C     ON SORT SI LAMATR
  369. C
  370.       IF(LMATR)GOTO 9998
  371. C
  372. C     DESACTIVATIONS
  373. C
  374.       SEGDES,MELEME,MPOVAL,MSOUPO
  375.       SEGSUP,MLENTI
  376.       SEGDES,MCHPOI
  377. C
  378. C     ON APPELLE RESO
  379. C
  380.       CALL RESOU
  381. C
  382. C     ... ET ON REGARDE LE RESULTAT
  383. C
  384.       CALL LIROBJ('CHPOINT',MCHPOI,1,IRETOU)
  385.       IF(IERR.NE.0)RETURN
  386. C
  387.       SEGACT,MCHPOI
  388.       MSOUPO=IPCHP(1)
  389.       SEGACT,MSOUPO
  390.       MELEME=IGEOC
  391.       MPOVAL=IPOVAL
  392.       SEGACT,MELEME,MPOVAL
  393. C
  394.       SEGINI,MLENTI,MLENT1
  395.       DO IE1=1,JG
  396.         LECT(IE1)=NUM(1,IE1)
  397.         MLENT1.LECT(IE1)=IE1
  398.       ENDDO
  399.       CALL GENOR2(LECT(1),MLENT1.LECT(1),JG)
  400. C
  401.       SEGINI,MLREEL
  402.       DO IE1=1,JG
  403.         PROG(IE1)=VPOCHA(MLENT1.LECT(IE1),1)
  404.       ENDDO
  405.       CALL VECTAB(MLREEL,JG,MTABLE)
  406. C
  407.       CALL ECROBJ('TABLE   ',MTABLE)
  408. C
  409. C     DESTRUCTION
  410. C
  411.       SEGSUP,MLENTI,MLENT1,MLREEL
  412. C
  413.       SEGSUP,MELEME,MPOVAL,MSOUPO,MCHPOI
  414. C
  415.       MCHPOI=MCHPO1
  416.       SEGACT,MCHPOI
  417.       MSOUPO=IPCHP(1)
  418.       SEGACT,MSOUPO
  419.       MELEME=IGEOC
  420.       MPOVAL=IPOVAL
  421.       SEGSUP,MELEME,MPOVAL,MSOUPO,MCHPOI
  422. C
  423.       SEGACT,MRIGID
  424.       MELEME=IRIGEL(1,1)
  425.       DESCR=IRIGEL(3,1)
  426.       xMATRI=IRIGEL(4,1)
  427. *      SEGACT,IMATRI
  428. *      XMATRI=IMATTT(1)
  429.       SEGSUP,MELEME,DESCR,XMATRI,MRIGID
  430. C
  431.       NBPTS=NBPTT
  432.       SEGADJ,MCOORD
  433. C
  434.       RETURN
  435. C
  436. C     ERREURS
  437. C
  438.  9998 SEGDES,MELEME,MPOVAL,MSOUPO
  439.       SEGSUP,MLENTI
  440.  9999 SEGDES,MCHPOI
  441.       RETURN
  442.       END
  443.  
  444.  
  445.  
  446.  
  447.  
  448.  
  449.  
  450.  
  451.  
  452.  
  453.  
  454.  
  455.  
  456.  
  457.  

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