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  26/02/23    21:15:05     12480          
  2.  
  3. *----------------------------------------------------------------------*
  4. *   CALCUL DE LA RIGIDITE POUR LA FORMULATION 'HHO'                    *
  5. *----------------------------------------------------------------------*
  6. *   ENTREES :                                                          *
  7. *   ---------                                                          *
  8. *                                                                      *
  9. *   SORTIES :                                                          *
  10. *   ---------                                                          *
  11. *----------------------------------------------------------------------*
  12.  
  13.       SUBROUTINE HHORIG (modHHO, IPRIGI,ISOU, IPDSCR,
  14.      &                   MATE,IVAMAT,NVMAT, IVACAR,NVCAR, iret)
  15.  
  16.       IMPLICIT INTEGER(I-N)
  17.       IMPLICIT REAL*8(A-H,O-Z)
  18.  
  19. -INC PPARAM
  20. -INC CCOPTIO
  21. -INC CCREEL
  22.  
  23. -INC CCHHOPA
  24. -INC CCHHOPR
  25.  
  26. -INC SMCOORD
  27. -INC SMCHAML
  28. -INC SMELEME
  29. c!!-INC SMINTE
  30. -INC SMLENTI
  31. -INC SMLREEL
  32.       POINTEUR mlrmsh.mlreel, mlrmbh.mlreel
  33. -INC SMMODEL
  34. -INC SMRIGID
  35.  
  36. -INC TMPTVAL
  37.  
  38.       SEGMENT MWKMAT
  39.         REAL*8 VALMAT(NTMAT), HOELAS(lhook,lhook)
  40. c**        REAL*8 XLOC(3,3),XGLOB(3,3),TXR(IDIM,IDIM)
  41.       ENDSEGMENT
  42.  
  43.       LOGICAL LDPGE
  44.  
  45.       iret = 0
  46.  
  47.       IPMATR = 0
  48.  
  49.       imodel = modHHO
  50. c*      segact,imodel <- actif en entree/sortie
  51.  
  52. C- Premieres verifications :
  53.       CALL HHONOB(modHHO, nobHHO, iret)
  54.       IF (nobHHO.LE.0) THEN
  55.         write(ioimp,*) 'HHORIG: data not consistent !'
  56.         iret = 5
  57.         call erreur(5)
  58.         return
  59.       END IF
  60.  
  61. C- Recuperation des donnees du modele en entree
  62.       meleme = imodel.IMAMOD
  63.       NBNO   = meleme.NUM(/1)
  64.       NBELEM = meleme.NUM(/2)
  65.  
  66. c*      MELE = imodel.NEFMOD
  67.       MFR  = imodel.infele(13)
  68.       CALL INFDPG(MFR,IFOUR,LDPGE,NDPGE)
  69.       IF (LDPGE) THEN
  70.         write(ioimp,*) 'HHORIG: PLAN GENE not implemented !'
  71.         iret = 21
  72.         RETURN
  73.       END IF
  74.  
  75. C---- !!
  76.       mlenti = imodel.IVAMOD(nobHHO+1)
  77. c*      segact,mlenti
  78.  
  79. c*      n_o_face = mlenti.lect(2)
  80. c*      n_d_face = mlenti.lect(3)
  81. c*      n_o_cell = mlenti.lect(4)
  82. c*      n_d_cell = mlenti.lect(5)
  83.       nb_faces = mlenti.lect(7)
  84.       NBPGAU   = mlenti.lect(8)
  85.       LRE      = mlenti.lect(11)
  86. c-dbg      write(ioimp,*) 'HHORIG:',ISOU,':',nb_faces,LRE,NBPGAU
  87.  
  88.       IF (nb_faces.NE.NBNO) THEN
  89.         write(ioimp,*) '(WARNING)HHORIG: Bizarre nb_faces'
  90.       END IF
  91.       IF (NBPGAU .NE. imodel.INFELE(6)) then
  92.         write(ioimp,*) '(WARNING)HHORIG: Bizarre nb.p.gau'
  93.       END IF
  94.       IF (LRE .NE. imodel.INFELE(9)) then
  95.         write(ioimp,*) '(WARNING)HHORIG: Bizarre LRE'
  96.       END IF
  97.  
  98. C- Verification des proprietes materielles : Champ CONSTANT par ELEMENT !
  99.       mptval = IVAMAT
  100.       NTMAT = mptval.IVAL(/1)
  101.       NUMAT = NTMAT - 4
  102.       lhook = 9
  103.  
  104. C- Tableau de stockage des proprietes materielles
  105.       MWKMAT = 0
  106.       SEGINI,MWKMAT
  107.       IWKMAT = MWKMAT
  108.  
  109.       isz = 0
  110.       iunif = 0
  111.       DO is = 1, NTMAT
  112.         IF (is.LE.(NVMAT-4) .OR. is.GT.NVMAT) THEN
  113.           melval = mptval.IVAL(is)
  114.           IGMN = melval.VELCHE(/1)
  115.           IEMN = melval.VELCHE(/2)
  116.           IF (IGMN.NE.1) isz = isz + 1
  117.           IF (IGMN.EQ.1 .AND. IEMN.EQ.1) iunif = iunif + 1
  118.           VALMAT(is) = melval.VELCHE(1,1)
  119.         END IF
  120.       END DO
  121.       IF (isz.NE.0) THEN
  122.         write(ioimp,*) 'HHORIG: Material properties must be constant ',
  123.      &                 'by element (for the moment) !'
  124.         iret = 21
  125.         RETURN
  126.       END IF
  127. * Le champ de proprietes est completement uniforme :
  128.       IF (iunif.NE.NUMAT) iunif = 0
  129.  
  130. C- Materiau isotrope et proprietes uniformes : -> calcul de HOELAS
  131.       IF (MATE.EQ.1 .AND. iunif.NE.0) THEN
  132.         CALL DOHMAS(VALMAT,'ISOTROPE',IFOUR,lhook,1,HOELAS,iretc)
  133.         r_z = HOELAS(4,4) * 2.D0
  134.         DO i = 4, lhook
  135.           HOELAS(i,i) = r_z
  136.         END DO
  137.       END IF
  138. C - Materiaux orthotrope et anisotrope
  139.       IF (MATE.EQ.2.OR.MATE.EQ.3) THEN
  140. C** AFAIRE !!!        
  141.       END IF
  142.  
  143. C- Les donnees liees a la formulation HHO :
  144.       IVCSTA = IVAL(NVMAT-3)
  145.       melval = IVCSTA
  146.       IGCS = melval.VELCHE(/1)
  147.       IECS = melval.VELCHE(/2)
  148. c-dbg      write(ioimp,*) 'IVCSTA',melval,igcs,iecs,tyval(nvmat-3)
  149.       IF (IGCS.NE.1) THEN
  150.         write(ioimp,*) 'HHORIG : Stabilization coefficient must be ',
  151.      &                 'constant by element'
  152.         iret = 21
  153.         RETURN
  154.       END IF
  155.  
  156.       IVMSTA = IVAL(NVMAT)
  157.       melval = IVMSTA
  158.       IGMS = melval.IELCHE(/1)
  159.       IEMS = melval.IELCHE(/2)
  160. c-dbg      write(ioimp,*) 'IVMSTA',melval,igms,iems,tyval(nvmat)
  161.       IF (IGMS.NE.1) THEN
  162.         write(ioimp,*) 'HHORIG : Stabilization matrix must be ',
  163.      &                 'constant by element'
  164.         iret = 21
  165.         RETURN
  166.       END IF
  167.       mlrmsh = melval.IELCHE(1,1)
  168. c*      segact,mlrmsh
  169. c*      write(ioimp,*) 'HHORIG MSTAB SIZE:',mlrmsh.prog(/1), LRE,LRE*LRE,
  170. c*     &               mlenti.lect(16)
  171.       IF ((mlrmsh.prog(/1).NE.(LRE*LRE)) .OR. 
  172.      &    (mlrmsh.prog(/1).NE.mlenti.lect(16))) THEN
  173.         write(ioimp,*) 'HHORIG : Stabilization matrix size incorrect'
  174.         iret = 21
  175.         RETURN
  176.       END IF
  177.  
  178.       IVPIHO = IVAL(NVMAT-2)
  179.       melval = IVPIHO
  180.       IGPI = melval.VELCHE(/1)
  181.       IEPI = melval.VELCHE(/2)
  182. c-dbg      write(ioimp,*) 'IVPIHO',melval,igpi,iepi,tyval(nvmat-2)
  183.  
  184.       IVMBHO = IVAL(NVMAT-1)
  185.       melval = IVMBHO
  186.       IGMB = melval.IELCHE(/1)
  187.       IEMB = melval.IELCHE(/2)
  188. c-dbg      write(ioimp,*) 'IVMBHO',melval,igmb,iemb,tyval(nvmat-1)
  189.       mlrmbh = melval.IELCHE(1,1)
  190. c*      segact,mlrmbh
  191. c*      write(ioimp,*) 'HHORIG MSTAB SIZE:',mlrmsh.prog(/1), LRE,LRE*LRE,
  192. c*     &               mlenti.lect(16)
  193.       IF ((mlrmbh.prog(/1).NE.(9*LRE)) .OR. 
  194.      &    (mlrmbh.prog(/1).NE.mlenti.lect(14))) THEN
  195.         write(ioimp,*) 'HHORIG : BHHO matrix size incorrect'
  196.         iret = 21
  197.         RETURN
  198.       END IF
  199.  
  200. C Le maillage support pour la rigidite est constitue du point support 
  201. C de la cellule et des points supports de chaque face de la cellule 
  202. C A FAIRE : ajouter le point support des DPGE !
  203.       NBNN = 1 + NBNO
  204.       IF (LDPGE) NBNN = NBNN + 1
  205.       NBREF  = 0
  206.       NBSOUS = 0
  207.       SEGINI,ipt1
  208.       ipt1.ITYPEL = 28
  209.       IF (LDPGE) THEN
  210.         DO i = 1, NBELEM
  211.           ipt1.NUM(NBNN,i) = IIPDPG
  212.         END DO
  213.       END IF
  214.       mlent1 = imodel.IVAMOD(nobHHO+4)
  215. c*      segact,mlent1
  216.       nbel1 = mlent1.lect(/1) / 2
  217.       if (nbel1.ne.NBELEM) write(ioimp,*) 'Bizarre hhorig nbel1'
  218.       ipt2 = MPCHHO
  219.       SEGACT,ipt2
  220.       DO i = 1, nbel1
  221.         je = mlent1.lect(2*i-1)
  222.         ip = mlent1.lect(2*i)
  223.         if (ip.le.0) write(ioimp,*) 'HHORIG: Bizarre...',i,je,ip
  224.         jp = ip + NBCHHO(je-1)
  225.         ipt1.num(1,i) = ipt2.num(1,jp)
  226.       END DO
  227.       mlent1 = imodel.IVAMOD(nobHHO+3)
  228. c*      segact,mlent1
  229.       nbfac1 = mlent1.lect(/1) / 2
  230. c-dbg      write(ioimp,*) 'HHORIG:',nbelem,nbno,nbfac1,nbnn
  231.       ipt2 = MPFHHO
  232.       SEGACT,ipt2
  233.       DO jg1 = 1, nbfac1
  234.         i_z = (jg1 - 1) / nbno
  235.         ie1 = i_z + 1
  236.         ip1 = jg1 - (i_z * nbno) + 1
  237.         je = mlent1.lect(2*jg1-1)
  238.         ip = ABS(mlent1.lect(2*jg1))
  239.         if (ip.eq.0) write(ioimp,*) 'HHORIG: Bizarre...',i,je,ip
  240.         jp = ip + NBFHHO(je-1)
  241.         ipt1.num(ip1,ie1) = ipt2.num(1,jp)
  242. c-dbg      write(ioimp,*) 'LO',jg1,ie1,ip1,'--',jp,je
  243. c-dbg      write(ioimp,*) '  ',ipt2.num(1,jp)
  244.       END DO
  245. c*      SEGDES,ipt1
  246.  
  247. C    INITIALISATION DU SEGMENT  XMATRI
  248.  
  249.       NLIGRP = LRE
  250.       NLIGRD = LRE
  251.       NELRIG = NBELEM
  252.       SEGINI,XMATRI
  253.       XMATRI.SYMRE = 0
  254.  
  255. C-------------------------
  256. C Boucle sur les cellules du modele
  257. C-------------------------
  258.       DO IEL = 1,NBELEM
  259.  
  260. C- Recuperation des proprietes constantes par element :
  261.         IF (iunif.EQ.0) THEN
  262.           DO is = 1, NVMAT
  263.             IF (is.LE.(NVMAT-4) .OR. is.GT.NVMAT) THEN
  264.               melval = mptval.IVAL(is)
  265.               VALMAT(is) = melval.VELCHE(1,IEL)
  266.             END IF
  267.           END DO
  268.         END IF
  269. C- Materiau isotrope
  270.         IF (MATE.EQ.1 .AND. iunif.EQ.0) THEN
  271.           CALL DOHMAS(VALMAT,'ISOTROPE',IFOUR,lhook,1,HOELAS,iretc)
  272.           r_z = HOELAS(4,4) * 2.D0
  273.           DO i = 4, lhook
  274.             HOELAS(i,i) = r_z
  275.           END DO
  276.         ENDIF
  277. C - Calcul des axes locaux dans les cas orthotrope et anisotrope
  278.         IF (MATE.EQ.2.OR.MATE.EQ.3) THEN
  279.         ENDIF
  280.  
  281. C- Initialisation de la matrice de rigidite avec la matrice de stabilisation
  282.         melval = IVCSTA
  283.         CSTAB  = melval.VELCHE(1,MIN(IEL,IECS))
  284.         melval = IVMSTA
  285.         mlrmsh = melval.IELCHE(1,MIN(IEL,IEMS))
  286. c*        segact,mlrmsh
  287. c* !! matrice SHHO stockee colonne par colonne : LRE*LRE (symetrique ?)
  288.         DO j = 1, LRE
  289.           jc = LRE * (j-1)
  290.           DO i = 1, LRE
  291.             RE(i,j,IEL) = CSTAB * mlrmsh.prog(jc+i)
  292.           END DO
  293.         END DO
  294. c*        segdes,mlrmsh
  295.  
  296.         JEPI = MIN(IEL,IEPI)
  297.         JEMB = MIN(IEL,IEMB)
  298. C--  --  --  --  --  --  --  --  --
  299. C - Boucle sur les points de Gauss
  300. C--  --  --  --  --  --  --  --  --
  301.         DO IGAU = 1, NBPGAU
  302.  
  303. C-AFAIRE C -- Recuperation des proprietes materielles (IGAU)
  304. C-AFAIRE           DO is = 1, NVMAT
  305. C-AFAIRE             IF (is.LE.(NVMAT-4) .OR. is.GT.NVMAT) THEN
  306. C-AFAIRE               melval = mptval.IVAL(is)
  307. C-AFAIRE               JEMN = MIN(IEL ,melval.VELCHE(/2))
  308. C-AFAIRE               JGMN = MIN(IGAU,melval.VELCHE(/1))
  309. C-AFAIRE               VALMAT(is) = melval.VELCHE(JGMN,JEMN)
  310. C-AFAIRE             END IF
  311. C-AFAIRE           END DO
  312. C-AFAIREC- Materiau isotrope
  313. C-AFAIRE           IF (MATE.EQ.1) THEN
  314. C-AFAIRE             CALL DOHMAS(VALMAT,'ISOTROPE',IFOUR,lhook,1,HOELAS,iretc)
  315. C-AFAIRE             r_z = HOELAS(4,4) * 2.D0
  316. C-AFAIRE             DO i = 4, lhook
  317. C-AFAIRE               HOELAS(i,i) = r_z
  318. C-AFAIRE             END DO
  319. C-AFAIRE           END IF
  320.  
  321.           melval = IVPIHO
  322.           JGPI = MIN(IGAU,IGPI)
  323.           poig = melval.VELCHE(JGPI,JEPI)
  324.  
  325.           melval = IVMBHO
  326.           JGMB = MIN(IGAU,IGMB)
  327.           mlrmbh = melval.IELCHE(JGMB,JEMB)
  328. c*          segact,mlrmbh
  329. c* !! matrice BHHO stockee colonne par colonne : lhook*LRE
  330.           CALL BDBST(mlrmbh.prog(1),poig,HOELAS,LRE,lhook,RE(1,1,IEL))
  331.  
  332. C--  --  --  --  --  --  --  --  --
  333.         END DO
  334. C--  --  --  --  --  --  --  --  --
  335. C-------------------------
  336.       END DO
  337. C-------------------------
  338.  
  339.       SEGDES,xMATRI
  340.       IPMATR = xMATRI
  341.  
  342.       mrigid = IPRIGI
  343. c*      segact,mrigid*mod
  344.  
  345.       mrigid.IRIGEL(1,ISOU) = ipt1
  346.       mrigid.IRIGEL(2,ISOU) = 0
  347.       mrigid.IRIGEL(3,ISOU) = IPDSCR
  348.       mrigid.IRIGEL(4,ISOU) = IPMATR
  349.       mrigid.IRIGEL(5,ISOU) = NIFOUR
  350.       mrigid.IRIGEL(6,ISOU) = 0
  351.       mrigid.IRIGEL(7,ISOU) = 0
  352.       mrigid.COERIG(ISOU) = 1.D0
  353.  
  354.       IF (MATE.EQ.2.OR.MATE.EQ.3) THEN
  355.       ENDIF
  356.       SEGSUP,MWKMAT
  357.  
  358. *      RETURN
  359.       END
  360.  
  361.  
  362.  

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