Télécharger la plaquette Cast3M


fr en

Télécharger frvisq.eso

Retour à la liste

Numérotation des lignes : 

frvisq

  1. C FRVISQ    SOURCE    MB234859  25/09/08    21:15:29     12358          
  2.  
  3.       SUBROUTINE FRVISQ(IPMODL,JPMAIL,IPCHE1, IPRIG)
  4. C
  5. C***********************************************************************
  6. C                                                                      *
  7. C  Routine principale appelée par AMOR                                 *
  8. C                                                                      *
  9. C  Calcule la matrice d'amortissement associée à la frontière du       *
  10. C  maillage dans plusieurs cas :                                       *
  11. C                                                                      *
  12. C   FORMULATION MECANIQUE                                              *
  13. C   +++++++++++++++++++++                                              *
  14. C                                                                      *
  15. C     * cas des massifs, dont l'enveloppe est constituée de SEG2 ou    *
  16. C       SEG3 (cas 2D), FAC3, FAC4, FAC6, ou FAC8 (cas 3D)              *
  17. C                                                                      *
  18. C   FORMULATION LIQUIDE                                                *
  19. C   +++++++++++++++++++                                                *
  20. C                                                                      *
  21. C     * cas des éléments dont l'enveloppe est constituée d'éléments    *
  22. C       à 2 (cas 2D), 3 ou 4 noeuds (cas 3D).                          *
  23. C______________________________________________________________________*
  24. C                                                                      *
  25. C  Entrées :                                                           *
  26. C  --------                                                            *
  27. C                                                                      *
  28. C     IPMODL : pointeur sur le modèle, objet MMODEL                    *
  29. C     JPMAIL : pointeur sur le maillage de la frontière, objet MELEME  *
  30. C     IPCHE1 : pointeur sur le champ par éléments de caractéristiques  *
  31. C              matériau, objet MCHAML                                  *
  32. C  Sorties :                                                           *
  33. C  --------                                                            *
  34. C     IPRIG : pointeur sur la matrice d'amortissement construite,      *
  35. C             objet MRIGID (=0 en cas d'erreur)                        *
  36. C                                                                      *
  37. C***********************************************************************
  38. C
  39.       IMPLICIT INTEGER(I-N)
  40.       IMPLICIT REAL*8(A-H,O-Z)
  41.  
  42. -INC PPARAM
  43. -INC CCOPTIO
  44. -INC CCHAMP
  45.  
  46. -INC SMCOORD
  47. -INC SMELEME
  48. -INC SMMODEL
  49. -INC SMCHAML
  50. -INC SMRIGID
  51. -INC SMINTE
  52.  
  53. -INC TMPTVAL
  54.  
  55.       INTEGER oooval
  56. C
  57.       SEGMENT NOTYPE
  58.         CHARACTER*16 TYPE(NBTYPE)
  59.       ENDSEGMENT
  60. C
  61.       CHARACTER*(NCONCH) CONM
  62.  
  63. C Support du champ de caracteristiques
  64.       PARAMETER ( IPLAZ=3 )
  65.  
  66.       PARAMETER (NINF=3)
  67.       INTEGER INFOS(NINF)
  68.  
  69.       IPRIG = 0
  70.  
  71.       IF (IFOUR.EQ.-3) THEN
  72.         CALL ERREUR(21)
  73.         RETURN
  74.       ENDIF
  75.  
  76. c_______________________________________________________________________
  77. c
  78. c     activation du modele
  79. c_______________________________________________________________________
  80. C
  81.       MMODEL=IPMODL
  82.       SEGACT MMODEL
  83.       NSOUS=KMODEL(/1)
  84.  
  85. C______________________________________________________________________C
  86. C                                                                      C
  87. C  CREATION DE L'OBJET MATRICE DE RIGIDITE                             C
  88. C______________________________________________________________________C
  89. C                                                                      C
  90.       NRIGEL=0
  91.       SEGINI,MRIGID
  92.       MTYMAT='RIGIDITE'
  93.       IFORIG=IFOUR
  94.       ICHOLE=0
  95.       IMGEO1=0
  96.       IMGEO2=0
  97.       ISUPEQ=0
  98.       JRCOND=0
  99.       JRDEPP=0
  100.       JRDEPD=0
  101.  
  102.       NFOR = 0
  103.       NMAT = 0
  104.       MN3  = 1
  105.       NOBMOD = 0
  106.       SEGINI,IMODE1
  107.       IMODE1.CMATEE = 'ISOTROPE'
  108.       MMODE1.KMODEL(1) = IMODE1
  109. C______________________________________________________________________C
  110. C                                                                      C
  111. C  BOUCLE SUR LES SOUS ZONES                                           C
  112. C______________________________________________________________________C
  113. C                                                                      C
  114.       DO 100 ISOUS = 1, NSOUS
  115. C
  116. C  on récupère l'information générale
  117. C
  118.         IMODEL = KMODEL(ISOUS)
  119.         SEGACT,IMODEL
  120.  
  121. C- Initialisations
  122.         IENVEL = 0
  123.         IPOGEO = 0
  124.  
  125.         IPT1 = IMAMOD
  126.         CONM = CONMOD
  127.         MELM = NEFMOD
  128. C
  129. C  création du tableau info
  130. C
  131.         iret = 1
  132.         CALL IDENT(IPT1,CONM,IPCHE1,0,INFOS,iret)
  133.         IF (iret.EQ.0) GOTO 1099
  134. C
  135. C  Determination de l'enveloppe du maillage massif du sous-modele
  136. C
  137.         CALL ECROBJ('MAILLAGE',IPT1)
  138.         IF (IDIM.EQ.3) THEN
  139.           CALL ENVELO
  140.         ELSE IF (IDIM.EQ.2) THEN
  141.           CALL PRCONT
  142.         ELSE
  143.           CALL ERREUR(5)
  144.         ENDIF
  145.         IF (IERR.NE.0) GOTO 1099
  146.         CALL LIROBJ('MAILLAGE',IENVEL,1,iret)
  147.         IF (IERR.NE.0) GOTO 1099
  148. C
  149. C  Elements de l'enveloppe IENVEL dans le maillage frontiere JPMAIL
  150. C
  151.         iret = 0
  152.         CALL INTERB(IENVEL,JPMAIL,iret,IPOGEO)
  153.         IF (iret.GT.0) GOTO 1099
  154.  
  155.         IPT3 = IPOGEO
  156.         SEGACT,IPT3
  157.         NBSOU3 = IPT3.LISOUS(/1)
  158.         IPT2 = IPT3
  159. C
  160.         NFOR = FORMOD(/2)
  161.         NMAT = MATMOD(/2)
  162.         MN3  = INFMOD(/1)
  163.         SEGADJ,IMODE1
  164.         IMODE1.CONMOD = CONM
  165.         DO i = 1, NFOR
  166.           IMODE1.FORMOD(i) = FORMOD(i)
  167.         ENDDO
  168.         DO i = 1, NMAT
  169.           IMODE1.MATMOD(i) = MATMOD(i)
  170.         ENDDO
  171.         DO i = 1, MN3
  172.           IMODE1.INFMOD(i) = INFMOD(i)
  173.         ENDDO
  174. C
  175. C  boucle sur les sous-zones de l'enveloppe
  176. C
  177.         DO 110 IB = 1, MAX(1,NBSOU3)
  178.  
  179. C-- Initialisations :
  180.           MOFORC = 0
  181.           MODEPL = 0
  182.           IPMINT = 0
  183.           MOMATR = 0
  184.           MOCARA = 0
  185.           MOTYPM = 0
  186.           MOTYPC = 0
  187.           ISUPM  = 0
  188.           ISUPC  = 0
  189.           IDESCR = 0
  190.  
  191. C-- Informations sur la (sous-zone de) l'enveloppe
  192.           IF (NBSOU3.NE.0) THEN
  193.             IPT2 = IPT3.LISOUS(IB)
  194.             SEGACT,IPT2
  195.           ENDIF
  196.           NBNOE2 = IPT2.NUM(/1)
  197.           NBELE2 = IPT2.NUM(/2)
  198.           LETYP  = IPT2.ITYPEL
  199. C-- Petit test sur le type
  200.           IF (LETYP.EQ.1) THEN
  201.             CALL ERREUR(16)
  202.             GOTO 1199
  203.           ENDIF
  204.           IPOGEO = IPT2
  205. C
  206. C-- On détermine la formulation associée à l'objet géométrique
  207. C-- elementaire de surface
  208.           CALL TYPFAC(MELM,NBNOE2,MELE)
  209. C
  210. C-- ERREUR : impossible d'utiliser FROABS pour les éléments
  211. C--          de formulation MELM
  212.           IF (MELE.EQ.0) THEN
  213.             MOTERR(1:8) = NOMTP(MELM)
  214.             CALL ERREUR(193)
  215.             GOTO 1199
  216.           ENDIF
  217.  
  218. C-- Information sur l'élément fini
  219. C
  220. C   Mise a jour de IMODE1 avec les donnees necessaires de IMODEL
  221. C
  222.           IMODE1.IMAMOD=IPOGEO
  223.           IMODE1.NEFMOD=MELE
  224.  
  225.           CALL PRQUOI(IMODE1)
  226.           IF (IERR.NE.0) GOTO 1199
  227. C
  228.           MFR   = IMODE1.INFELE(13)
  229.           LRE   = IMODE1.INFELE(9)
  230.           LW    = IMODE1.INFELE(7)
  231.           NDDL  = IMODE1.INFELE(15)
  232.           IPPORE = 0
  233.           IPMINT = IMODE1.INFELE(11)
  234.  
  235. C-- Recherche des inconnues primales et duales (DEPL-FORC)
  236.           CALL IDPRIM(IMODEL,MFR,MODEPL,NDEPL,ndum)
  237.           CALL IDDUAL(IMODEL,MFR,MOFORC,NFORC,ndum)
  238.  
  239.           IF (NDEPL.EQ.0 .OR. NFORC.EQ.0 .OR. NDEPL.NE.NFORC) THEN
  240.             CALL ERREUR(5)
  241.             GOTO 1199
  242.           ENDIF
  243.  
  244. C-- Remplissage du segment DESCRipteur
  245.           NLIGRP = LRE
  246.           NLIGRD = LRE
  247.           SEGINI,DESCR
  248.  
  249.           NCOMP = NDEPL
  250.           NBNNS = NBNOE2
  251.           IDDL=1
  252.           DO 1004 INOEUD=1,NBNNS
  253.             DO 1005 ICOMP=1,NCOMP
  254.               NOMID = MODEPL
  255.               LISINC(IDDL)=LESOBL(ICOMP)
  256.               NOMID = MOFORC
  257.               LISDUA(IDDL)=LESOBL(ICOMP)
  258.               NOELEP(IDDL)=INOEUD
  259.               NOELED(IDDL)=INOEUD
  260.               IDDL=IDDL+1
  261.  1005       CONTINUE
  262.  1004     CONTINUE
  263.  
  264.           IDESCR = DESCR
  265.  
  266. C-- Recuperation des noms de composantes MATERIAU
  267.           nbrobl = 0
  268.           nbrfac = 0
  269.           nomid  = 0
  270.           notype = 0
  271.  
  272. C  rho, E, nu pour les massifs
  273.           IF (MFR.EQ.1) THEN
  274.             nbrobl = 3
  275.             SEGINI,nomid
  276.             lesobl(1) = 'RHO '
  277.             lesobl(2) = 'YOUN'
  278.             lesobl(3) = 'NU  '
  279.  
  280.             nbtype = 1
  281.             SEGINI,notype
  282.             type(1) = 'REAL*8'
  283. C
  284. C  rho, cson, rhoref, cref, rlcar pour les liquides
  285.           ELSE IF (MFR.EQ.11.OR.MFR.EQ.41) THEN
  286.             nbrobl = 5
  287.             SEGINI,nomid
  288.             lesobl(1) = 'RHO '
  289.             lesobl(2) = 'CSON'
  290.             lesobl(3) = 'RORF'
  291.             lesobl(4) = 'CREF'
  292.             lesobl(5) = 'LCAR'
  293.  
  294.             nbtype = 1
  295.             SEGINI,notype
  296.             type(1) = 'REAL*8'
  297.           ENDIF
  298.  
  299.           MOMATR = nomid
  300.           MOTYPM = notype
  301.           NMATR = nbrobl
  302.           NMATF = nbrfac
  303.           NMATT = NMATR+NMATF
  304.  
  305. C--- Verification du support des composantes recherchées
  306.           IF (MOMATR.NE.0) THEN
  307.             CALL QUESUQ(IMODEL,IPCHE1,3,0,MOMATR,IPLAZ,ISUPM,iret)
  308.             IF (ISUPM.GT.1) GOTO 1199
  309.           ENDIF
  310. C
  311. C-- Recuperation des noms de composantes CARACTERISTIQUES
  312.           nbrobl = 0
  313.           nbrfac = 0
  314.           nomid  = 0
  315.           notype = 0
  316.  
  317. C Epaisseur du massif en contraintes planes
  318.           IF (MFR.EQ.1 .AND. IFOUR.EQ.-2) THEN
  319.             nbrfac = 1
  320.             SEGINI,nomid
  321.             lesfac(1) = 'DIM3'
  322.             nbtype = 1
  323.             SEGINI,notype
  324.             type(1) = 'REAL*8'
  325.           ENDIF
  326.  
  327.           MOCARA = nomid
  328.           MOTYPC = notype
  329.           NCARA = nbrobl
  330.           NCARF = nbrfac
  331.           NCARR = NCARA+NCARF
  332.  
  333. C--- Verification du support des composantes recherchées
  334.           IF (MOCARA.NE.0) THEN
  335.             CALL QUESUQ(IMODEL,IPCHE1,3,0,MOCARA,IPLAZ,ISUPC,iret)
  336.             IF (ISUPC.GT.1) GOTO 1199
  337.           ENDIF
  338.  
  339. C-- Segment d'integration MINTE
  340.           MINTE = IPMINT
  341.           SEGACT,MINTE
  342.           NBPGAU = POIGAU(/1)
  343.  
  344. C- Partionnement si necessaire de la matrice d'amortissement
  345. C- determinant ainsi le nombre d'objets elementaires de MRIGID
  346.           LTRK = oooval(1,4)
  347.           IF (LTRK.EQ.0) LTRK = oooval(1,1)
  348.       LTRK=MAX(LTRK,2**24)
  349. * Ajout a la taille en mots de la matrice des infos du segment
  350.           LSEG = LRE*LRE*NBELE2 + 16
  351.           NBLPRT = (LSEG-1)/LTRK + 1
  352.           NBLMAX = (NBELE2-1)/NBLPRT + 1
  353.           NBLPRT = (NBELE2-1)/NBLMAX + 1
  354. c*        write(ioimp,*) ' frvisq : nblprt nblmax = ',nblprt,nblmax,nbele2
  355. C*OF : Pour l'instant pas de partition pour FRVISQ
  356.           NBLPRT = 1
  357.  
  358. C-- Ajout de la matrice d'AMORTISSEMENT a la matrice globale
  359.           NRIGE0 = IRIGEL(/2)
  360.           NRIGEL = NRIGE0 + NBLPRT
  361.           SEGADJ,MRIGID
  362.  
  363.           descr  = IDESCR
  364.           meleme = IPOGEO
  365.           nbnn   = NBNOE2
  366.           nbelem = NBELE2
  367.           nbsous = 0
  368.           nbref  = 0
  369.  
  370.           DO 120 irige = 1, NBLPRT
  371.  
  372. C-- Mettre ici la partition du maillage IPOGEO
  373.             ipmail = meleme
  374.             ipdesc = descr
  375.  
  376. C- Initialisation de la matrice de rigidite elementaire (xmatri)
  377.             NELRIG = nbelem
  378.             SEGINI,xmatri
  379.             ipmatr = xmatri
  380.  
  381. C- Recuperation des valeurs des proprietes materiau et geometriques
  382. c* Note : les proprietes sont les valeurs au support des EF massifs
  383. c*        et non celles au niveau de l'enveloppe surfacique !
  384. c*        Cela ne marche que si les proprietes sont constantes. Dans
  385. c*        les autres cas, le resultat est... Pour eviter cela, on met
  386. c*        un test sur la constance du champ !
  387.             ivamat = 0
  388.             ivacar = 0
  389.             IF (MOMATR.NE.0) THEN
  390.               CALL KOMCHA(IPCHE1,IPT1,CONM,MOMATR,MOTYPM,1,
  391. c*              CALL KOMCHA(IPCHE1,ipmail,CONM,MOMATR,MOTYPM,1,
  392.      &                    INFOS,3, ivamat)
  393.               IF (IERR.NE.0) GOTO 1199
  394.               mptval = ivamat
  395.               do i = 1, NMATT
  396.                 if (ival(i).ne.0) then
  397.                   melval = IVAL(i)
  398.                   if (velche(/1).ne.1 .and. velche(/2).ne.1) then
  399.                     write(ioimp,*) 'Champ MATERIAU non constant'
  400.                     call erreur(21)
  401.                     goto 1199
  402.                   endif
  403.                 endif
  404.               enddo
  405.               IF (ISUPM.EQ.1) THEN
  406.                 CALL VALCHE(ivamat,NMATT,IPMINT,IPPORE,MOMATR,MELE)
  407.                 IF (IERR.NE.0) THEN
  408.                   ISUPM = 0
  409.                   GOTO 1199
  410.                 ENDIF
  411.               ENDIF
  412.             ENDIF
  413.             IF (MOCARA.NE.0) THEN
  414.               CALL KOMCHA(IPCHE1,IPT1,CONM,MOCARA,MOTYPC,1,
  415. c*              CALL KOMCHA(IPCHE1,ipmail,CONM,MOCARA,MOTYPC,1,
  416.      &                    INFOS,3, ivacar)
  417.               IF (IERR.NE.0) GOTO 1199
  418.               mptval = ivacar
  419.               do i = 1, NCARR
  420.                 if (ival(i).ne.0) then
  421.                   melval = IVAL(i)
  422.                   if (velche(/1).ne.1 .and. velche(/2).ne.1) then
  423.                     write(ioimp,*) 'Champ MATERIAU non constant'
  424.                     call erreur(21)
  425.                     goto 1199
  426.                   endif
  427.                 endif
  428.               enddo
  429.               IF (ISUPC.EQ.1) THEN
  430.                 CALL VALCHE(ivacar,NCARR,IPMINT,IPPORE,MOCARA,MELE)
  431.                 IF (IERR.NE.0)THEN
  432.                   ISUPC = 0
  433.                   GOTO 1199
  434.                 ENDIF
  435.               ENDIF
  436.             ENDIF
  437.  
  438. C  distinction des cas 2D et 3D
  439. C______________________________________________________________________C
  440. C                                                                      C
  441. C            CAS DES ELEMENTS MASSIFS BIDIMENSIONNELS                  C
  442. C                  FACES ASSOCIEES SEG2 OU SEG3                        C
  443. C______________________________________________________________________C
  444. C                                                                      C
  445.          IF (MELE.EQ.2.OR.MELE.EQ.3) THEN
  446. C
  447.             CALL FROA2D(ipmail,ipmatr,IPMINT,ivamat,ivacar,
  448.      1                  MELE,MFR,LRE,NDDL)
  449. C
  450. C______________________________________________________________________C
  451. C                                                                      C
  452. C               CAS DES ELEMENTS LIQUIDES 2D OU 3D                     C
  453. C               FACES ASSOCIEES LSE2, LTR3 OU LQU4                     C
  454. C______________________________________________________________________C
  455. C                                                                      C
  456. C
  457.          ELSE IF(MELE.EQ.97.OR.MELE.EQ.35.OR.MELE.EQ.36) THEN
  458. C
  459.             CALL LFROA(ipmail,ipmatr,IPMINT,ivamat,ivacar,
  460.      1                  MELE,MFR,LRE,NDDL)
  461. C
  462. C______________________________________________________________________C
  463. C                                                                      C
  464. C            CAS DES ELEMENTS MASSIFS TRIDIMENSIONNELS                 C
  465. C             FACES ASSOCIEES FAC3,FAC4,FAC6 OU FAC8                   C
  466. C______________________________________________________________________C
  467. C
  468.          ELSE IF(MELE.EQ.31.OR.MELE.EQ.32.OR.MELE.EQ.33.
  469.      1                                   OR.MELE.EQ.34)THEN
  470. C
  471.             CALL FROA3D(ipmail,ipmatr,IPMINT,ivamat,ivacar,
  472.      1                  MELE,MFR,LRE,NDDL)
  473. C
  474. C  erreur, l'élément n'est pas encore implémenté
  475. C
  476.          ELSE
  477. C
  478.             MOTERR(1:4)=NOMTP(MELE)
  479.             MOTERR(5:12)='FRVISQ'
  480.             CALL ERREUR (86)
  481.          ENDIF
  482. C
  483.          IF (ISUPM.EQ.1 .OR. nblprt.GT.1) THEN
  484.            CALL DTMVAL(ivamat,3)
  485.          ELSE
  486.            CALL DTMVAL(ivamat,1)
  487.          ENDIF
  488.          IF (ISUPC.EQ.1 .OR. nblprt.GT.1) THEN
  489.            CALL DTMVAL(ivacar,3)
  490.          ELSE
  491.            CALL DTMVAL(ivacar,1)
  492.          ENDIF
  493.  
  494.           xmatri = ipmatr
  495.           IF (NBLPRT.GT.1) THEN
  496.             meleme = ipmail
  497.           ENDIF
  498.  
  499. C- Sortie prematuree en cas d'erreur
  500.           IF (IERR.NE.0) GOTO 1199
  501.  
  502. C- Stockage de la matrice
  503.           jrige = NRIGE0 + irige
  504.           COERIG(jrige)   = 1.
  505.           IRIGEL(1,jrige) = ipmail
  506.           IRIGEL(2,jrige) = 0
  507.           IRIGEL(3,jrige) = ipdesc
  508.           IRIGEL(4,jrige) = ipmatr
  509.           IRIGEL(5,jrige) = NIFOUR
  510.           IRIGEL(6,jrige) = 0
  511.           IRIGEL(7,jrige) = 0
  512. * matrice non symetrique (forces sur pi seulement
  513. * qui dependent de p)
  514.           IF (MFR.EQ.11.OR.MFR.EQ.41) THEN
  515.             IRIGEL(7,jrige) = 2
  516.             xmatri.symre=2
  517.            ENDIF
  518.           SEGDES,xmatri
  519.           IRIGEL(8,jrige) = 0
  520.  
  521.  120      CONTINUE
  522. C- Fin de la boucle de partition maillage/rigidite
  523.  
  524.  1199     CONTINUE
  525.           IF (MOMATR.NE.0) THEN
  526.             nomid = MOMATR
  527.             SEGSUP,nomid
  528.             notype = MOTYPM
  529.             SEGSUP,notype
  530.           ENDIF
  531.           IF (MOCARA.NE.0) THEN
  532.             nomid = MOCARA
  533.             SEGSUP,nomid
  534.             notype = MOTYPC
  535.             SEGSUP,notype
  536.           ENDIF
  537. C
  538.  
  539. C- Sortie prematuree en cas d'erreur
  540.           IF (IERR.NE.0) GOTO 1098
  541.  
  542.  110    CONTINUE
  543. C- Fin de la boucle sur (les sous-zones de) l'enveloppe
  544. C
  545.  1098   CONTINUE
  546.  1099   CONTINUE
  547. C- Sortie prematuree en cas d'erreur
  548.         IF (IERR.NE.0) GOTO 999
  549. C
  550.  100  CONTINUE
  551.       SEGSUP,IMODE1
  552. C- Fin de la boucle sur les modeles elementaires
  553. C
  554.       NRIGE0 = IRIGEL(/2)
  555.       IF (NRIGE0.EQ.0) THEN
  556.         CALL ERREUR(902)
  557.       ENDIF
  558.  
  559.  999  CONTINUE
  560.       IF (IERR.EQ.0) THEN
  561.         IPRIG = MRIGID
  562.         SEGDES,MRIGID
  563.       ELSE
  564.         IPRIG = 0
  565.         SEGSUP,MRIGID
  566.       ENDIF
  567.  
  568. c      RETURN
  569.       END
  570.  
  571.  
  572.  
  573.  

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