Télécharger la plaquette Cast3M


fr en

Télécharger hhorig.eso

Retour à la liste

Numérotation des lignes : 

hhorig

  1. C HHORIG    SOURCE    OF166741  24/05/06    21:15:14     11082          
  2. C HHORIG    SOURCE    FANDEUR                                           
  3.  
  4. *----------------------------------------------------------------------*
  5. *   CALCUL DE LA RIGIDITE POUR LA FORMULATION 'HHO'                    *
  6. *----------------------------------------------------------------------*
  7. *   ENTREES :                                                          *
  8. *   ________                                                           *
  9. *        MATE     Numero du materiau                                   *
  10. *        IVAMAT   Pointeur sur un segment MPTVAL pour le materiau ou   *
  11. *        NVMAT    Nombre de composante de materiau                     *
  12. *                                                                      *
  13. *   SORTIES :                                                          *
  14. *   ________                                                           *
  15. *        IPMATR   pointeur sur la rigidite de la sous-zone             *
  16. *----------------------------------------------------------------------*
  17.  
  18.       SUBROUTINE HHORIG (modHHO, IPRIGI,ISOU,
  19.      &                   MATE,IVAMAT,NVMAT, IVACAR,NVCAR, iret)
  20.  
  21.       IMPLICIT INTEGER(I-N)
  22.       IMPLICIT REAL*8(A-H,O-Z)
  23.  
  24. -INC PPARAM
  25. -INC CCOPTIO
  26. -INC CCREEL
  27.  
  28. -INC CCHHOPA
  29. -INC CCHHOPR
  30.  
  31. -INC SMCOORD
  32. -INC SMCHAML
  33. -INC SMELEME
  34. c!!-INC SMINTE
  35. -INC SMLENTI
  36. -INC SMLREEL
  37.       POINTEUR mlrmsh.mlreel, mlrmbh.mlreel
  38. -INC SMMODEL
  39.       POINTEUR nomidu.nomid, nomidf.nomid
  40. -INC SMRIGID
  41.  
  42.       SEGMENT MWKMAT
  43.         REAL*8 VALMAT(NTMAT), HOELAS(lhook,lhook)
  44. c**        REAL*8 XLOC(3,3),XGLOB(3,3),TXR(IDIM,IDIM)
  45.       ENDSEGMENT
  46.  
  47.       SEGMENT MPTVAL
  48.         INTEGER IPOS(NS),NSOF(NS)
  49.         INTEGER IVAL(NCOSOU)
  50.         CHARACTER*16 TYVAL(NCOSOU)
  51.       ENDSEGMENT
  52.  
  53.       LOGICAL LDPGE
  54.  
  55.       iret = 0
  56.  
  57.       MRIGID = IPRIGI
  58. c*      segact,MRIGID*mod
  59.  
  60.       IPMATR = 0
  61.       IPDSCR = 0
  62.  
  63.       imodel = modHHO
  64. c*      segact,imodel <- actif en entree/sortie
  65.  
  66. C- Premieres verifications :
  67.       CALL HHONOB(modHHO, nobHHO, iret)
  68.       IF (nobHHO.LE.0) THEN
  69.         write(ioimp,*) 'HHORIG: data not consistent !'
  70.         call erreur(5)
  71.         return
  72.       END IF
  73.  
  74. C- Recuperation des donnees de infell en entree
  75.       meleme = imodel.IMAMOD
  76.       MELE  = imodel.NEFMOD
  77.       MFR   = imodel.infele(13)
  78.       CALL INFDPG(MFR,IFOUR,LDPGE,NDPGE)
  79.       IF (LDPGE) THEN
  80.         write(ioimp,*) 'HHORIG: PLAN GENE not implemented !'
  81.         iret = 21
  82.         RETURN
  83.       END IF
  84.  
  85. C---- !!
  86.       mlenti = imodel.IVAMOD(nobHHO+1)
  87. c*      segact,mlenti
  88.  
  89. c*      n_o_face = mlenti.lect(2)
  90.       n_d_face = mlenti.lect(3)
  91. c*      n_o_cell = mlenti.lect(4)
  92.       n_d_cell = mlenti.lect(5)
  93.       nb_faces = mlenti.lect(7)
  94.       NBPGAU   = mlenti.lect(8)
  95.       idifo    = mlenti.lect(9)
  96.       LRE      = mlenti.lect(11)
  97.       NDEPF    = idifo * n_d_face
  98.       NDEPC    = idifo * n_d_cell
  99.       NBDDL    = LRE
  100. c-dbg      write(ioimp,*) 'HHORIG:',ISOU,':',nb_faces,NBDDL,NBPGAU
  101.  
  102.       IF (nb_faces.NE.meleme.NUM(/1)) THEN
  103.         write(ioimp,*) 'HHORIG: Bizarre nb_faces'
  104.       END IF
  105.       IF (LRE .NE. imodel.INFELE(9)) then
  106.         write(ioimp,*) 'HHORIG: Bizarre LRE'
  107.       END IF
  108.       IF (NBPGAU .NE. imodel.INFELE(6)) then
  109.         write(ioimp,*) 'HHORIG: Bizarre nb.p.gau'
  110.       END IF
  111.  
  112. C- Verification des proprietes materielles : Champ CONSTANT par ELEMENT !
  113.       mptval = IVAMAT
  114.       NTMAT = mptval.IVAL(/1)
  115.       NUMAT = NTMAT - 4
  116.       lhook = 9
  117.  
  118. C- Tableau de stockage des proprietes materielles
  119.       MWKMAT = 0
  120.       SEGINI,MWKMAT
  121.       IWKMAT = MWKMAT
  122.  
  123.       isz = 0
  124.       iunif = 0
  125.       DO is = 1, NTMAT
  126.         IF (is.LE.(NVMAT-4) .OR. is.GT.NVMAT) THEN
  127.           melval = mptval.IVAL(is)
  128.           IGMN = melval.VELCHE(/1)
  129.           IEMN = melval.VELCHE(/2)
  130.           IF (IGMN.NE.1) isz = isz + 1
  131.           IF (IGMN.EQ.1 .AND. IEMN.EQ.1) iunif = iunif + 1
  132.           VALMAT(is) = melval.VELCHE(1,1)
  133.         END IF
  134.       END DO
  135.       IF (isz.NE.0) THEN
  136.         write(ioimp,*) 'HHORIG: Material properties must be constant ',
  137.      &                 'by element (for the moment) !'
  138.         iret = 21
  139.         RETURN
  140.       END IF
  141. * Le champ de proprietes est completement uniforme :
  142.       IF (iunif.NE.NUMAT) iunif = 0
  143.  
  144. C- Materiau isotrope et proprietes uniformes : -> calcul de HOELAS
  145.       IF (MATE.EQ.1 .AND. iunif.NE.0) THEN
  146.         CALL DOHMAS(VALMAT,'ISOTROPE',IFOUR,lhook,1,HOELAS,iretc)
  147.         r_z = HOELAS(4,4) * 2.D0
  148.         DO i = 4, lhook
  149.           HOELAS(i,i) = r_z
  150.         END DO
  151.       END IF
  152. C - Materiaux orthotrope et anisotrope
  153.       IF (MATE.EQ.2.OR.MATE.EQ.3) THEN
  154. C** AFAIRE !!!        
  155.       END IF
  156.  
  157. C- Les donnees liees a la formulation HHO :
  158.       IVCSTA = IVAL(NVMAT-3)
  159.       melval = IVCSTA
  160.       IGCS = melval.VELCHE(/1)
  161.       IECS = melval.VELCHE(/2)
  162. c-dbg      write(ioimp,*) 'IVCSTA',melval,igcs,iecs,tyval(nvmat-3)
  163.       IF (IGCS.NE.1) THEN
  164.         write(ioimp,*) 'HHORIG : Stabilization coefficient must be ',
  165.      &                 'constant by element'
  166.         iret = 21
  167.         RETURN
  168.       END IF
  169.  
  170.       IVMSTA = IVAL(NVMAT)
  171.       melval = IVMSTA
  172.       IGMS = melval.IELCHE(/1)
  173.       IEMS = melval.IELCHE(/2)
  174. c-dbg      write(ioimp,*) 'IVMSTA',melval,igms,iems,tyval(nvmat)
  175.       IF (IGMS.NE.1) THEN
  176.         write(ioimp,*) 'HHORIG : Stabilization matrix must be ',
  177.      &                 'constant by element'
  178.         iret = 21
  179.         RETURN
  180.       END IF
  181.       mlrmsh = melval.IELCHE(1,1)
  182. c*      segact,mlrmsh
  183. c*      write(ioimp,*) 'HHORIG MSTAB SIZE:',mlrmsh.prog(/1), LRE,LRE*LRE,
  184. c*     &               mlenti.lect(16)
  185.       IF ((mlrmsh.prog(/1).NE.(LRE*LRE)) .OR. 
  186.      &    (mlrmsh.prog(/1).NE.mlenti.lect(16))) THEN
  187.         write(ioimp,*) 'HHORIG : Stabilization matrix size incorrect'
  188.         iret = 21
  189.         RETURN
  190.       END IF
  191.  
  192.       IVPIHO = IVAL(NVMAT-2)
  193.       melval = IVPIHO
  194.       IGPI = melval.VELCHE(/1)
  195.       IEPI = melval.VELCHE(/2)
  196. c-dbg      write(ioimp,*) 'IVPIHO',melval,igpi,iepi,tyval(nvmat-2)
  197.  
  198.       IVMBHO = IVAL(NVMAT-1)
  199.       melval = IVMBHO
  200.       IGMB = melval.IELCHE(/1)
  201.       IEMB = melval.IELCHE(/2)
  202. c-dbg      write(ioimp,*) 'IVMBHO',melval,igmb,iemb,tyval(nvmat-1)
  203.       mlrmbh = melval.IELCHE(1,1)
  204. c*      segact,mlrmbh
  205. c*      write(ioimp,*) 'HHORIG MSTAB SIZE:',mlrmsh.prog(/1), LRE,LRE*LRE,
  206. c*     &               mlenti.lect(16)
  207.       IF ((mlrmbh.prog(/1).NE.(9*LRE)) .OR. 
  208.      &    (mlrmbh.prog(/1).NE.mlenti.lect(14))) THEN
  209.         write(ioimp,*) 'HHORIG : BHHO matrix size incorrect'
  210.         iret = 21
  211.         RETURN
  212.       END IF
  213.  
  214. C Le maillage support pour la rigidite est constitue du point support 
  215. C de la cellule et des points supports de chaque face de la cellule 
  216. C A FAIRE : ajouter le point support des DPGE !
  217.       NBNO   = meleme.NUM(/1)
  218.       NBNN   = 1 + NBNO
  219.       IF (LDPGE) NBNN = NBNN + 1
  220.       NBELEM = meleme.NUM(/2)
  221.       NBREF  = 0
  222.       NBSOUS = 0
  223.       SEGINI,ipt1
  224.       ipt1.ITYPEL = 28
  225.       IF (LDPGE) THEN
  226.         DO i = 1, NBELEM
  227.           ipt1.NUM(NBNN,i) = IIPDPG
  228.         END DO
  229.       END IF
  230.       mlent1 = imodel.IVAMOD(nobHHO+4)
  231. c*      segact,mlent1
  232.       nbel1 = mlent1.lect(/1) / 2
  233.       if (nbel1.ne.NBELEM) write(ioimp,*) 'Bizarre hhorig nbel1'
  234.       ipt2 = MPCHHO
  235.       SEGACT,ipt2
  236.       DO i = 1, nbel1
  237.         je = mlent1.lect(2*i-1)
  238.         ip = mlent1.lect(2*i)
  239.         if (ip.le.0) write(ioimp,*) 'HHORIG: Bizarre...',i,je,ip
  240.         jp = ip + NBCHHO(je-1)
  241.         ipt1.num(1,i) = ipt2.num(1,jp)
  242.       END DO
  243.       mlent1 = imodel.IVAMOD(nobHHO+3)
  244. c*      segact,mlent1
  245.       nbfac1 = mlent1.lect(/1) / 2
  246. c-dbg      write(ioimp,*) 'HHORIG:',nbelem,nbno,nbfac1,nbnn
  247.       ipt2 = MPFHHO
  248.       SEGACT,ipt2
  249.       DO jg1 = 1, nbfac1
  250.         i_z = (jg1 - 1) / nbno
  251.         ie1 = i_z + 1
  252.         ip1 = jg1 - (i_z * nbno) + 1
  253.         je = mlent1.lect(2*jg1-1)
  254.         ip = ABS(mlent1.lect(2*jg1))
  255.         if (ip.eq.0) write(ioimp,*) 'HHORIG: Bizarre...',i,je,ip
  256.         jp = ip + NBFHHO(je-1)
  257.         ipt1.num(ip1,ie1) = ipt2.num(1,jp)
  258. c-dbg      write(ioimp,*) 'LO',jg1,ie1,ip1,'--',jp,je
  259. c-dbg      write(ioimp,*) '  ',ipt2.num(1,jp)
  260.       END DO
  261. c*      SEGDES,ipt1
  262.  
  263. C  REMPLISSAGE DU SEGMENT DESCRIPTEUR
  264. C
  265. c!         NSTRS = INFELE(16)
  266. c!         LW    = INFELE(7)
  267. c!         NDDL  = INFELE(15)
  268. c!         LHOOK = INFELE(10)
  269.  
  270. C        RECHERCHE DES NOMS D'INCONNUES ET DES DUAUX
  271. C
  272.       nomidu = imodel.lnomid(1)
  273.       ndepl  = nomidu.lesobl(/2)
  274. c-      ndum   = nomidu.lesfac(/2)
  275.  
  276.       nomidf = imodel.lnomid(2)
  277.       nforc  = nomidf.lesobl(/2)
  278. c-      ndum   = nomidf.lesfac(/2)
  279.  
  280.       IF (NDEPL.EQ.0.OR.NFORC.EQ.0.OR.NDEPL.NE.NFORC) THEN
  281.         write(ioimp,*) 'HHORIG: data inconsistent PRIMAL/DUAL (1)'
  282.         iret = 5
  283.       END IF
  284.       IF (NDEPL.NE.(NDEPC+NDEPF)) THEN
  285.         write(ioimp,*) 'HHORIG: data inconsistent PRIMAL/DUAL (2)'
  286.         iret = 717
  287.       END IF
  288.       IF (NBDDL.NE.(NDEPC+(nb_faces*NDEPF))) THEN
  289.         write(ioimp,*) 'HHORIG: data inconsistent PRIMAL/DUAL (3)'
  290.         iret = 717
  291.       END IF
  292.       IF (iret.NE.0) RETURN
  293.  
  294.       NLIGRP = LRE
  295.       NLIGRD = LRE
  296.       SEGINI,DESCR
  297.  
  298.       iddl = 0
  299. c  Noeud support de la cellule :
  300.       INOE = 1
  301.       DO ic = 1, NDEPC
  302.         iddl = iddl + 1
  303.         LISINC(iddl) = nomidu.LESOBL(ic)
  304.         LISDUA(iddl) = nomidf.LESOBL(ic)
  305.         NOELEP(iddl) = INOE
  306.         NOELED(iddl) = INOE
  307.       END DO
  308. c  Noeuds supports des faces de la cellule :
  309.       DO in = 1, NBNO
  310.         INOE = in + 1
  311.         DO ic = 1, NDEPF
  312.           iddl = iddl + 1
  313.           LISINC(iddl) = nomidu.LESOBL(NDEPC+ic)
  314.           LISDUA(iddl) = nomidf.LESOBL(NDEPC+ic)
  315.           NOELEP(iddl) = INOE
  316.           NOELED(iddl) = INOE
  317.         END DO
  318.       END DO
  319. C CAS DE LA DEFORMATION PLANE GENERALISEE
  320. c!!      IF (LDPGE) THEN
  321. C!!        INOE = NBNN
  322. c!!        DO ic = (NDPGE-1),0,-1
  323. c!!          iddl = iddl + 1
  324. c!!          LISINC(iddl) = nomidu.LESOBL(NDEPL-ic)
  325. c!!          LISDUA(iddl) = nomidf.LESOBL(NFORC-ic)
  326. c!!          NOELEP(iddl) = INOE
  327. c!!          NOELED(iddl) = INOE
  328. c!!        END DO
  329. c!!      END IF
  330.       IF (iddl .NE. NBDDL) THEN
  331.         write(ioimp,*) 'HHORIG: NBDDL inconsistent',iddl,nbddl
  332.         iret = 717
  333.         RETURN
  334.       END IF
  335.  
  336.       SEGDES,DESCR
  337.       IPDSCR = DESCR
  338.  
  339. C    INITIALISATION DU SEGMENT  XMATRI
  340. C
  341.       NLIGRP = LRE
  342.       NLIGRD = LRE
  343.       NELRIG = NBELEM
  344.       SEGINI,XMATRI
  345.       XMATRI.SYMRE = 0
  346.  
  347. C-------------------------
  348. C Boucle sur les cellules du modele
  349. C-------------------------
  350.       DO IEL = 1,NBELEM
  351.  
  352. C- Recuperation des proprietes constantes par element :
  353.         IF (iunif.EQ.0) THEN
  354.           DO is = 1, NVMAT
  355.             IF (is.LE.(NVMAT-4) .OR. is.GT.NVMAT) THEN
  356.               melval = mptval.IVAL(is)
  357.               VALMAT(is) = melval.VELCHE(1,IEL)
  358.             END IF
  359.           END DO
  360.         END IF
  361. C- Materiau isotrope
  362.         IF (MATE.EQ.1 .AND. iunif.EQ.0) THEN
  363.           CALL DOHMAS(VALMAT,'ISOTROPE',IFOUR,lhook,1,HOELAS,iretc)
  364.           r_z = HOELAS(4,4) * 2.D0
  365.           DO i = 4, lhook
  366.             HOELAS(i,i) = r_z
  367.           END DO
  368.         ENDIF
  369. C - Calcul des axes locaux dans les cas orthotrope et anisotrope
  370.         IF (MATE.EQ.2.OR.MATE.EQ.3) THEN
  371.         ENDIF
  372.  
  373. C- Initialisation de la matrice de rigidite avec la matrice de stabilisation
  374.         melval = IVCSTA
  375.         CSTAB  = melval.VELCHE(1,MIN(IEL,IECS))
  376.         melval = IVMSTA
  377.         mlrmsh = melval.IELCHE(1,MIN(IEL,IEMS))
  378. c*        segact,mlrmsh
  379. c* !! matrice SHHO stockee colonne par colonne : LRE*LRE (symetrique ?)
  380.         DO j = 1, LRE
  381.           jc = LRE * (j-1)
  382.           DO i = 1, LRE
  383.             RE(i,j,IEL) = CSTAB * mlrmsh.prog(jc+i)
  384.           END DO
  385.         END DO
  386. c*        segdes,mlrmsh
  387.  
  388.         JEPI = MIN(IEL,IEPI)
  389.         JEMB = MIN(IEL,IEMB)
  390. C--  --  --  --  --  --  --  --  --
  391. C - Boucle sur les points de Gauss
  392. C--  --  --  --  --  --  --  --  --
  393.         DO IGAU = 1, NBPGAU
  394.  
  395. C-AFAIRE C -- Recuperation des proprietes materielles (IGAU)
  396. C-AFAIRE           DO is = 1, NVMAT
  397. C-AFAIRE             IF (is.LE.(NVMAT-4) .OR. is.GT.NVMAT) THEN
  398. C-AFAIRE               melval = mptval.IVAL(is)
  399. C-AFAIRE               JEMN = MIN(IEL ,melval.VELCHE(/2))
  400. C-AFAIRE               JGMN = MIN(IGAU,melval.VELCHE(/1))
  401. C-AFAIRE               VALMAT(is) = melval.VELCHE(JGMN,JEMN)
  402. C-AFAIRE             END IF
  403. C-AFAIRE           END DO
  404. C-AFAIREC- Materiau isotrope
  405. C-AFAIRE           IF (MATE.EQ.1) THEN
  406. C-AFAIRE             CALL DOHMAS(VALMAT,'ISOTROPE',IFOUR,lhook,1,HOELAS,iretc)
  407. C-AFAIRE             r_z = HOELAS(4,4) * 2.D0
  408. C-AFAIRE             DO i = 4, lhook
  409. C-AFAIRE               HOELAS(i,i) = r_z
  410. C-AFAIRE             END DO
  411. C-AFAIRE           END IF
  412.  
  413.           melval = IVPIHO
  414.           JGPI = MIN(IGAU,IGPI)
  415.           poig = melval.VELCHE(JGPI,JEPI)
  416.  
  417.           melval = IVMBHO
  418.           JGMB = MIN(IGAU,IGMB)
  419.           mlrmbh = melval.IELCHE(JGMB,JEMB)
  420. c*          segact,mlrmbh
  421. c* !! matrice BHHO stockee colonne par colonne : lhook*LRE
  422.           CALL BDBST(mlrmbh.prog(1),poig,HOELAS,LRE,lhook,RE(1,1,IEL))
  423.  
  424. C--  --  --  --  --  --  --  --  --
  425.         END DO
  426. C--  --  --  --  --  --  --  --  --
  427. C-------------------------
  428.       END DO
  429. C-------------------------
  430.  
  431.       SEGDES,XMATRI
  432.       IPMATR = XMATRI
  433.  
  434.       IRIGEL(1,ISOU) = ipt1
  435.       IRIGEL(2,ISOU) = 0
  436.       IRIGEL(3,ISOU) = IPDSCR
  437.       IRIGEL(4,ISOU) = IPMATR
  438.       IRIGEL(5,ISOU) = NIFOUR
  439.       IRIGEL(6,ISOU) = 0
  440.       IRIGEL(7,ISOU) = 0
  441.       COERIG(ISOU) = 1.D0
  442.  
  443.       IF (MATE.EQ.2.OR.MATE.EQ.3) THEN
  444.       ENDIF
  445.       SEGSUP,MWKMAT
  446.  
  447. *      RETURN
  448.       END
  449.  
  450.  
  451.  

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