Télécharger la plaquette Cast3M


fr en

Télécharger hhosig.eso

Retour à la liste

Numérotation des lignes : 

hhosig

  1. C HHOSIG    SOURCE    OF166741  25/02/21    21:17:27     12166          
  2.  
  3. *----------------------------------------------------------------------*
  4. *   CALCUL DES CONTRAINTES POUR LA FORMULATION 'HHO'                   *
  5. *----------------------------------------------------------------------*
  6. *   ENTREES :                                                          *
  7. *        MATE     Numero du materiau                                   *
  8. *        IVAMAT   Pointeur sur un segment MPTVAL pour le materiau      *
  9. *        NVMAT    Nombre de composantes obligatoires materiau          *
  10. *                                                                      *
  11. *   SORTIES :                                                          *
  12. *----------------------------------------------------------------------*
  13.  
  14.       SUBROUTINE HHOSIG(imoHHO, ichDEP,nmoDEP,
  15.      &                  IIPDPG,UZDPG,RYDPG,RXDPG,
  16.      &                  MATE,IVAMAT,NVMAT, IPMINT, NBPTEL,
  17.      &                  IVASTR,NCSTR, iret)
  18.  
  19.       IMPLICIT INTEGER(I-N)
  20.       IMPLICIT REAL*8(A-H,O-Z)
  21.  
  22. -INC PPARAM
  23. -INC CCOPTIO
  24. -INC CCREEL
  25.  
  26. -INC CCHHOPA
  27. -INC CCHHOPR
  28.  
  29. -INC SMCOORD
  30. -INC SMCHAML
  31. -INC SMELEME
  32. -INC SMINTE
  33. -INC SMLENTI
  34.       POINTEUR mlent4.mlenti
  35. -INC SMLMOTS
  36. -INC SMLREEL
  37.       POINTEUR mlrdef.mlreel, mlrdec.mlreel, mlrmbh.mlreel
  38. -INC SMMODEL
  39.  
  40. -INC TMPTVAL
  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 MWKHHO
  48.         INTEGER TABINT(NBINT)
  49.         REAL*8  TABFLO(NBFLO)
  50.       ENDSEGMENT
  51.  
  52.       DIMENSION UDPGE(3)
  53.  
  54.       iret = 0
  55.  
  56.       imodel = imoHHO
  57. c*      segact,imodel <- actif en entree/sortie
  58.  
  59. C- Premieres verifications :
  60.       CALL HHONOB(imoHHO, nobHHO, iret)
  61.       IF (nobHHO.LE.0) THEN
  62.         write(ioimp,*) 'HHOSIG: IMODEL incorrect (not HHO)'
  63.         iret = 5
  64.         return
  65.       END IF
  66.  
  67. C- Introduction du point autour duquel se fait le mouvement
  68. C  de la section en defo plane generalisee
  69. C  IIPDPG = numero du noeud/point support si defini pour le modele
  70. C  NDPGE > 0 si prise en compte du point support
  71.       IF (IIPDPG.GT.0) THEN
  72.         write(ioimp,*) 'HHOSIG: 2D PLAN GENE mode not implemented!'
  73.         iret = 5
  74.         return
  75.         IF (IFOUR.EQ.-3) THEN
  76.           NDPGE = 3
  77.           UDPGE(1) = UZDPG
  78.           UDPGE(2) = RYDPG
  79.           UDPGE(3) = RXDPG
  80. C          SEGACT,MCOORD
  81.           IREF = (IIPDPG-1)*(IDIM+1)
  82.           XDPGE = mcoord.XCOOR(IREF+1)
  83.           YDPGE = mcoord.XCOOR(IREF+2)
  84.         ELSE IF (IFOUR.EQ.11) THEN
  85.           NDPGE = 2
  86.           UDPGE(1) = UZDPG
  87.           UDPGE(2) = RXDPG
  88.           UDPGE(3) = XZero
  89.           XDPGE = XZero
  90.           YDPGE = XZero
  91.         ELSE IF (IFOUR.EQ. 7 .OR. IFOUR.EQ. 8 .OR. IFOUR.EQ. 9 .OR.
  92.      &           IFOUR.EQ.10 .OR. IFOUR.EQ.14) THEN
  93.           NDPGE = 1
  94.           UDPGE(1) = UZDPG
  95.           UDPGE(2) = XZero
  96.           UDPGE(3) = XZero
  97.           XDPGE = XZero
  98.           YDPGE = XZero
  99.         else
  100.           write(ioimp,*) 'HHOSIG: ERREUR NDPGE'
  101.           iret = 5
  102.           return
  103.         END IF
  104.       ELSE
  105.         NDPGE = 0
  106.         UDPGE(1) = XZero
  107.         UDPGE(2) = XZero
  108.         UDPGE(3) = XZero
  109.         XDPGE = XZero
  110.         YDPGE = XZero
  111.       END IF
  112.  
  113. C- Recuperation des donnees du modele en entree
  114. c*      MELE   = imodel.NEFMOD
  115. c*      MFR    = imodel.infele(13)
  116.       meleme = imodel.IMAMOD
  117.       NBNOE  = meleme.NUM(/1)
  118.       NBELT  = meleme.NUM(/2)
  119.  
  120.       mlenti = imodel.IVAMOD(nobHHO+1)
  121. c*      segact,mlenti
  122.       mlent3 = imodel.IVAMOD(nobHHO+3)
  123. c*      segact,mlent3
  124.       mlent4 = imodel.IVAMOD(nobHHO+4)
  125. c*      segact,mlent4
  126.  
  127. c*      n_o_face = mlenti.lect(2)
  128.       n_d_face = mlenti.lect(3)
  129. c*      n_o_cell = mlenti.lect(4)
  130.       n_d_cell = mlenti.lect(5)
  131.       nb_faces = mlenti.lect(7)
  132.       NBPGAU   = mlenti.lect(8)
  133.       idifo    = mlenti.lect(9)
  134.       NBDDL    = mlenti.lect(11)
  135.       NDEPF    = idifo * n_d_face
  136.       NDEPC    = idifo * n_d_cell
  137.       lhook = 9
  138.  
  139.       nbel3 = mlent3.lect(/1) / 2
  140.       nbfa3 = 2 * nb_faces
  141.       nbel4 = mlent4.lect(/1) / 2
  142.  
  143.       IF (mlenti.lect(6).NE.NBNOE) THEN
  144.         write(ioimp,*) 'HHOSIG: Bizarre nb_vertices'
  145.       END IF
  146.       IF (NBDDL .NE. imodel.INFELE(9)) then
  147.         write(ioimp,*) 'HHOSIG: Bizarre NBDDL'
  148.       END IF
  149. c NBPGAU =? (NBPTEL = imodel.INFELE(4))
  150.       IF (NBPGAU .NE. imodel.INFELE(4)) then
  151.         write(ioimp,*) 'HHOSRIG: Bizarre nb.p.gau(1)'
  152.       END IF
  153. c NBPGAU =? minte.POIGAU(/1)
  154.       minte = IPMINT
  155. c*      SEGACT minte <- actif en E/S
  156.       if (NBPGAU .NE. minte.POIGAU(/1)) then
  157.         write(ioimp,*) 'HHOSIG: Bizarre nb.p.gau (2)'
  158.       end if
  159. c-dbg      write(ioimp,*) 'HHOSIG:',nb_faces,NBDDL,NBPGAU
  160.       if (nbel3.NE.(NBELT*nb_faces)) then
  161.         write(ioimp,*) 'HHOSIG: Bizarre nbel3'
  162.       end if
  163.       if (nbel4.NE.NBELT) then
  164.         write(ioimp,*) 'HHOSIG: Bizarre nbel4'
  165.       end if
  166.  
  167.       ivid = 1
  168. C- Deplacements des faces et des cellules :
  169.       nomid = nmoDEP
  170. c*      segact,nomid
  171.       JGN = nomid.lesobl(/1)
  172. c-dbg      write(ioimp,*) 'HHOSIG=',NDEPC,NDEPF,lesobl(/2),lesfac(/2)
  173.       nfac = 0
  174. C Deplacements des cellules - Points supports des cellules
  175.       JGM = NDEPC + nfac
  176.       SEGINI,mlmots
  177.       DO i = 1, NDEPC
  178.         mlmots.MOTS(i)(1:JGN) = nomid.lesobl(i)(1:JGN)
  179.       END DO
  180. c*      DO i = 1, nfac
  181. c*        mlmots.MOTS(NDEP+i)(1:JGN) = nomid.lesfac(i)(1:JGN)
  182. c*      END DO
  183.       CALL EXTR23(ichDEP,mlmots,MPCHHO,mlrDEC,ivid)
  184.       SEGACT,mlrDEC
  185. c*      IF (IERR.NE.0) THEN
  186. c*        iret = 21
  187. c*        return
  188. c*      END IF
  189.  
  190. c-dbg      write(ioimp,*) 'mlmots DEPC ',(mots(i),i=1,mots(/2))
  191. c-dbg      write(ioimp,*) 'U.CELL',mlrdec.prog(/1),NCEHHO,NDEPC
  192. c-dbg      write(ioimp,*) '      ',(mlrdec.prog(i),i=1,mlrdec.prog(/1))
  193. c-dbg      write(ioimp,*)
  194.  
  195. C  Deplacements des faces - Points supports des faces
  196.       JGM = NDEPF + nfac
  197.       SEGADJ,mlmots
  198.       DO i = 1, NDEPF
  199.         mlmots.MOTS(i)(1:JGN) = nomid.lesobl(NDEPC+i)(1:JGN)
  200.       END DO
  201. c*      DO i = 1, nfac
  202. c*        mlmots.MOTS(NDEP+i)(1:JGN) = nomid.lesfac(i)(1:JGN)
  203. c*      END DO
  204.       CALL EXTR23(ichDEP,mlmots,MPFHHO,mlrDEF,ivid)
  205.       SEGACT,mlrDEF
  206. c*      IF (IERR.NE.0) THEN
  207. c*        iret = 21
  208. c*        return
  209. c*      END IF
  210.  
  211. c-dbg      write(ioimp,*) 'MLMOTS DEPF ',(mots(i),i=1,mots(/2))
  212. c-dbg      write(ioimp,*) 'U.FACE',mlrDEF.prog(/1),NFAHHO,NDEPF
  213. c-dbg      write(ioimp,*) '      ',(mlrDEF.prog(i),i=1,mlrDEF.prog(/1))
  214. c-dbg      write(ioimp,*)
  215.  
  216. C- Verification des proprietes materielles : Champ CONSTANT par ELEMENT !
  217.       mptval = IVAMAT
  218.       NTMAT = mptval.IVAL(/1)
  219.       NUMAT = NTMAT - 4
  220.  
  221. C- Tableau de stockage des proprietes materielles
  222.       MWKMAT = 0
  223.       SEGINI,MWKMAT
  224.       IWKMAT = MWKMAT
  225.  
  226.       isz = 0
  227.       iunif = 0
  228.       DO is = 1, NTMAT
  229.         IF (is.LE.(NVMAT-4) .OR. is.GT.NVMAT) THEN
  230.           melval = mptval.IVAL(is)
  231.           IGMN = melval.VELCHE(/1)
  232.           IEMN = melval.VELCHE(/2)
  233.           IF (IGMN.NE.1) isz = isz + 1
  234.           IF (IGMN.EQ.1 .AND. IEMN.EQ.1) iunif = iunif + 1
  235.           VALMAT(is) = melval.VELCHE(1,1)
  236.         END IF
  237.       END DO
  238.       IF (isz.NE.0) THEN
  239.         write(ioimp,*) 'HHOSIG: Material properties must be constant ',
  240.      &                 'by element (for the moment) !'
  241.         iret = 21
  242.         RETURN
  243.       END IF
  244. * Le champ de proprietes est completement uniforme :
  245.       IF (iunif.NE.NUMAT) iunif = 0
  246.  
  247. C- Materiau isotrope et proprietes uniformes : -> calcul de HOELAS
  248.       IF (MATE.EQ.1 .AND. iunif.NE.0) THEN
  249.         CALL DOHMAS(VALMAT,'ISOTROPE',IFOUR,lhook,1,HOELAS,iretc)
  250.         r_z = HOELAS(4,4) * 2.D0
  251.         DO i = 4, lhook
  252.           HOELAS(i,i) = r_z
  253.         END DO
  254.       END IF
  255. C - Materiaux orthotrope et anisotrope
  256.       IF (MATE.EQ.2.OR.MATE.EQ.3) THEN
  257. C** AFAIRE !!!        
  258.       END IF
  259.  
  260. C- Les donnees liees a la formulation HHO :
  261.       IVMBHO = mptval.IVAL(NVMAT-1)
  262.       melval = IVMBHO
  263.       IGMB = melval.IELCHE(/1)
  264.       IEMB = melval.IELCHE(/2)
  265. c-dbg      write(ioimp,*) 'IVMBHO',melval,igmb,iemb,tyval(nvmat-1)
  266.       mlrmbh = melval.IELCHE(1,1)
  267. c*      segact,mlrmbh
  268. c*      write(ioimp,*) 'HHOSIG MBHHO SIZE:',mlrmbh.prog(/1),
  269. c*     &                NBDDL,9*NBDDL,mlenti.lect(14)
  270.       IF ((mlrmbh.prog(/1).NE.(9*NBDDL)) .OR. 
  271.      &    (mlrmbh.prog(/1).NE.mlenti.lect(14))) THEN
  272.         write(ioimp,*) 'HHOSIG: BHHO matrix size incorrect'
  273.         iret = 21
  274.         RETURN
  275.       END IF
  276.  
  277. C- Indices et tableau de travail
  278.       ir_coo =          0
  279. c* si besoin des coordonnees      ir_eps = ir_coo + (IDIM*NBNOE)
  280.       ir_eps = ir_coo + 0
  281.       ir_sig = ir_eps + lhook
  282.       ir_uce = ir_sig + lhook
  283.       ir_ufa = ir_uce + NDEPC
  284.       ir_uge = ir_ufa + (NDEPF*nb_faces)
  285.       ir_fin = ir_uge + NDPGE
  286.  
  287.       NBINT = 1
  288.       NBFLO = ir_fin
  289.       SEGINI,MWKHHO
  290.  
  291.       IF (NDPGE.GT.0) THEN
  292.         DO ic = 1, NDPGE
  293.           TABFLO(ir_uge+ic) = UDPGE(ic)
  294.         END DO
  295.       END IF
  296. c* si besoin des coordonnees        SEGACT,MCOORD*NOMOD
  297.  
  298. C-------------------------
  299. C Boucle sur les cellules du modele
  300. C-------------------------
  301.       DO IEL = 1,NBELT
  302.  
  303. C- Recuperation des coordonnees des noeuds de l element IEL
  304. c*        CALL HHOCOO(meleme.num,NBNOE, mcoord.xcoor, IEL,
  305. c*    &               TABFLO(ir_coo+1), iret)
  306. c*        IF (iret.NE.0) RETURN
  307.  
  308. C- Valeurs des inconnues primales pour l'element IEL (cell+faces)
  309.         in1 = IEL * 2
  310.         je = mlent4.lect(in1-1)
  311.         ip = mlent4.lect(in1)
  312.       if (ip.le.0) write(ioimp,*) 'HHOSIG ICEL Bizarre...',iel,je,ip
  313.         jp = ip + NBCHHO(je-1)
  314.         DO ic = 1, NDEPC
  315.           jc = NCEHHO * (ic - 1)
  316.           TABFLO(ir_uce + ic) = mlrDEC.prog(jp + jc)
  317.         END DO
  318. c-dbg      write(ioimp,*) 'HHOSIG ICEL',iel,je,ip,jp,ir_uce
  319. c-dbg      write(ioimp,*) (TABFLO(ir_uce+ic),ic=1,ndepc)
  320.  
  321.         in1 = (IEL-1) * nbfa3
  322.         ir_kc = ir_ufa
  323.         DO j1 = 1, nb_faces
  324.           in2 = in1 + (2 * j1)
  325.           je = mlent3.lect(in2-1)
  326.           ip = ABS(mlent3.lect(in2))
  327.       if (ip.eq.0) write(ioimp,*) 'HHOSIG IFAE Bizarre...',iel,j1,je,ip
  328.           jp = ip + NBFHHO(je-1)
  329.           DO ic = 1, NDEPF
  330.             jc = NFAHHO * (ic - 1)
  331.             TABFLO(ir_kc + ic) = mlrDEF.prog(jp + jc)
  332.           END DO
  333. c-dbg      write(ioimp,*) 'HHOSIG IFAE',iel,j1,je,ip,jp,ir_kc
  334. c-dbg      write(ioimp,*) (TABFLO(ir_kc+ic),ic=1,ndepf)
  335.           ir_kc = ir_kc + NDEPF
  336.         END DO
  337.  
  338. c-dbg      write(ioimp,*) 'HHOSIG ...',ir_uce,ir_kc,NBDDL
  339.  
  340. C- Recuperation des proprietes constantes par element :
  341.         IF (iunif.EQ.0) THEN
  342.           DO is = 1, NVMAT
  343.             IF (is.LE.(NVMAT-4) .OR. is.GT.NVMAT) THEN
  344.               melval = mptval.IVAL(is)
  345.               VALMAT(is) = melval.VELCHE(1,IEL)
  346.             END IF
  347.           END DO
  348.         END IF
  349. C- Materiau isotrope
  350.         IF (MATE.EQ.1 .AND. iunif.EQ.0) THEN
  351.           CALL DOHMAS(VALMAT,'ISOTROPE',IFOUR,lhook,1,HOELAS,iretc)
  352.           r_z = HOELAS(4,4) * 2.D0
  353.           DO i = 4, lhook
  354.             HOELAS(i,i) = r_z
  355.           END DO
  356.         ENDIF
  357. C - Calcul des axes locaux dans les cas orthotrope et anisotrope
  358.         IF (MATE.EQ.2.OR.MATE.EQ.3) THEN
  359.         ENDIF
  360.  
  361.         JEMB = MIN(IEL,IEMB)
  362. C--  --  --  --  --  --  --  --  --
  363. C - Boucle sur les points de Gauss
  364. C--  --  --  --  --  --  --  --  --
  365.         DO IGAU = 1, NBPGAU
  366.  
  367. C-AFAIRE C -- Recuperation des proprietes materielles (IGAU)
  368. C-AFAIRE           DO is = 1, NVMAT
  369. C-AFAIRE             IF (is.LE.(NVMAT-4) .OR. is.GT.NVMAT) THEN
  370. C-AFAIRE               melval = mptval.IVAL(is)
  371. C-AFAIRE               JEMN = MIN(IEL ,melval.VELCHE(/2))
  372. C-AFAIRE               JGMN = MIN(IGAU,melval.VELCHE(/1))
  373. C-AFAIRE               VALMAT(is) = melval.VELCHE(JGMN,JEMN)
  374. C-AFAIRE             END IF
  375. C-AFAIRE           END DO
  376. C-AFAIREC- Materiau isotrope
  377. C-AFAIRE           IF (MATE.EQ.1) THEN
  378. C-AFAIRE             CALL DOHMAS(VALMAT,'ISOTROPE',IFOUR,lhook,1,HOELAS,iretc)
  379. C-AFAIRE             r_z = HOELAS(4,4) * 2.D0
  380. C-AFAIRE             DO i = 4, lhook
  381. C-AFAIRE               HOELAS(i,i) = r_z
  382. C-AFAIRE             END DO
  383. C-AFAIRE           END IF
  384.  
  385.           melval = IVMBHO
  386.           JGMB = MIN(IGAU,IGMB)
  387.           mlrmbh = melval.IELCHE(JGMB,JEMB)
  388. c*          segact,mlrmbh
  389. c* !! matrice BHHO stockee colonne par colonne : lhook*NBDDL
  390.           DO ic = 1, lhook
  391.             r_z = XZero
  392.             DO jc = 1, NBDDL
  393.               jnc = lhook * (jc-1)
  394.               r_z = r_z + mlrmbh.prog(ic + jnc) * TABFLO(ir_uce + jc)
  395.             END DO
  396.             TABFLO(ir_eps + ic) = r_z
  397.           END DO
  398.           DO ic = 1, lhook
  399.             r_z = XZero
  400.             DO jc = 1, lhook
  401.               r_z = r_z + HOELAS(ic,jc) * TABFLO(ir_eps + jc)
  402.             END DO
  403.             TABFLO(ir_sig + ic) = r_z
  404.           END DO
  405.  
  406. C -- Remplissage du segment contenant les "contraintes"
  407. c* A voir selon ordre - permutation :
  408.           mptval = IVASTR
  409.           DO ic = 1, NCSTR
  410.             melval = mptval.IVAL(ic)
  411.             melval.velche(IGAU,IEL) = TABFLO(ir_sig + ic)
  412.           END DO
  413.  
  414. C--  --  --  --  --  --  --  --  --
  415.         END DO
  416. C--  --  --  --  --  --  --  --  --
  417. C-------------------------
  418.       END DO
  419. C-------------------------
  420.  
  421.       IF (MATE.EQ.2.OR.MATE.EQ.3) THEN
  422.       ENDIF
  423.       SEGSUP,MWKMAT,MWKHHO
  424.       SEGSUP,mlrDEC,mlrDEF
  425.       SEGSUP,mlmots
  426.  
  427. c*      RETURN
  428.       END
  429.  
  430.  
  431.  

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