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

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