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    PV090527  26/04/28    21:15:35     12529          
  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. C______________________________________________________________________C
  109. C                                                                      C
  110. C  BOUCLE SUR LES SOUS ZONES                                           C
  111. C______________________________________________________________________C
  112. C                                                                      C
  113.       DO 100 ISOUS = 1, NSOUS
  114. C
  115. C  on récupère l'information générale
  116. C
  117.         IMODEL = KMODEL(ISOUS)
  118.         SEGACT,IMODEL
  119.  
  120. C- Initialisations
  121.         IENVEL = 0
  122.         IPOGEO = 0
  123.  
  124.         IPT1 = IMAMOD
  125.         CONM = CONMOD
  126.         MELM = NEFMOD
  127. C
  128. C  création du tableau info
  129. C
  130.         iret = 1
  131.         CALL IDENT(IPT1,CONM,IPCHE1,0,INFOS,iret)
  132.         IF (iret.EQ.0) GOTO 1099
  133. C
  134. C  Determination de l'enveloppe du maillage massif du sous-modele
  135. C
  136.         CALL ECROBJ('MAILLAGE',IPT1)
  137.         IF (IDIM.EQ.3) THEN
  138.           CALL ENVELO
  139.         ELSE IF (IDIM.EQ.2) THEN
  140.           CALL PRCONT
  141.         ELSE
  142.           CALL ERREUR(5)
  143.         ENDIF
  144.         IF (IERR.NE.0) GOTO 1099
  145.         CALL LIROBJ('MAILLAGE',IENVEL,1,iret)
  146.         IF (IERR.NE.0) GOTO 1099
  147. C
  148. C  Elements de l'enveloppe IENVEL dans le maillage frontiere JPMAIL
  149. C
  150.         iret = 0
  151.         CALL INTERB(IENVEL,JPMAIL,iret,IPOGEO)
  152.         IF (iret.GT.0) GOTO 1099
  153.  
  154.         IPT3 = IPOGEO
  155.         SEGACT,IPT3
  156.         NBSOU3 = IPT3.LISOUS(/1)
  157.         IPT2 = IPT3
  158. C
  159.         NFOR = FORMOD(/2)
  160.         NMAT = MATMOD(/2)
  161.         MN3  = INFMOD(/1)
  162.         SEGADJ,IMODE1
  163.         IMODE1.CONMOD = CONM
  164.         DO i = 1, NFOR
  165.           IMODE1.FORMOD(i) = FORMOD(i)
  166.         ENDDO
  167.         DO i = 1, NMAT
  168.           IMODE1.MATMOD(i) = MATMOD(i)
  169.         ENDDO
  170.         DO i = 1, MN3
  171.           IMODE1.INFMOD(i) = INFMOD(i)
  172.         ENDDO
  173. C
  174. C  boucle sur les sous-zones de l'enveloppe
  175. C
  176.         DO 110 IB = 1, MAX(1,NBSOU3)
  177.  
  178. C-- Initialisations :
  179.           MOFORC = 0
  180.           MODEPL = 0
  181.           IPMINT = 0
  182.           MOMATR = 0
  183.           MOCARA = 0
  184.           MOTYPM = 0
  185.           MOTYPC = 0
  186.           ISUPM  = 0
  187.           ISUPC  = 0
  188.           IDESCR = 0
  189.  
  190. C-- Informations sur la (sous-zone de) l'enveloppe
  191.           IF (NBSOU3.NE.0) THEN
  192.             IPT2 = IPT3.LISOUS(IB)
  193.             SEGACT,IPT2
  194.           ENDIF
  195.           NBNOE2 = IPT2.NUM(/1)
  196.           NBELE2 = IPT2.NUM(/2)
  197.           LETYP  = IPT2.ITYPEL
  198. C-- Petit test sur le type
  199.           IF (LETYP.EQ.1) THEN
  200.             CALL ERREUR(16)
  201.             GOTO 1199
  202.           ENDIF
  203.           IPOGEO = IPT2
  204. C
  205. C-- On détermine la formulation associée à l'objet géométrique
  206. C-- elementaire de surface
  207.           CALL TYPFAC(MELM,NBNOE2,MELE)
  208. C
  209. C-- ERREUR : impossible d'utiliser FROABS pour les éléments
  210. C--          de formulation MELM
  211.           IF (MELE.EQ.0) THEN
  212.             MOTERR(1:8) = NOMTP(MELM)
  213.             CALL ERREUR(193)
  214.             GOTO 1199
  215.           ENDIF
  216.  
  217. C-- Information sur l'élément fini
  218. C
  219. C   Mise a jour de IMODE1 avec les donnees necessaires de IMODEL
  220. C
  221.           IMODE1.IMAMOD=IPOGEO
  222.           IMODE1.NEFMOD=MELE
  223.  
  224.           CALL PRQUOI(IMODE1)
  225.           IF (IERR.NE.0) GOTO 1199
  226. C
  227.           MFR   = IMODE1.INFELE(13)
  228.           LRE   = IMODE1.INFELE(9)
  229.           LW    = IMODE1.INFELE(7)
  230.           NDDL  = IMODE1.INFELE(15)
  231.           IPPORE = 0
  232.           IPMINT = IMODE1.INFELE(11)
  233.  
  234. C-- Recherche des inconnues primales et duales (DEPL-FORC)
  235.           CALL IDPRIM(IMODEL,MFR,MODEPL,NDEPL,ndum)
  236.           CALL IDDUAL(IMODEL,MFR,MOFORC,NFORC,ndum)
  237.  
  238.           IF (NDEPL.EQ.0 .OR. NFORC.EQ.0 .OR. NDEPL.NE.NFORC) THEN
  239.             CALL ERREUR(5)
  240.             GOTO 1199
  241.           ENDIF
  242.  
  243. C-- Remplissage du segment DESCRipteur
  244.           NLIGRP = LRE
  245.           NLIGRD = LRE
  246.           SEGINI,DESCR
  247.  
  248.           NCOMP = NDEPL
  249.           NBNNS = NBNOE2
  250.           IDDL=1
  251.           DO 1004 INOEUD=1,NBNNS
  252.             DO 1005 ICOMP=1,NCOMP
  253.               NOMID = MODEPL
  254.               LISINC(IDDL)=LESOBL(ICOMP)
  255.               NOMID = MOFORC
  256.               LISDUA(IDDL)=LESOBL(ICOMP)
  257.               NOELEP(IDDL)=INOEUD
  258.               NOELED(IDDL)=INOEUD
  259.               IDDL=IDDL+1
  260.  1005       CONTINUE
  261.  1004     CONTINUE
  262.  
  263.           IDESCR = DESCR
  264.  
  265. C-- Recuperation des noms de composantes MATERIAU
  266.           nbrobl = 0
  267.           nbrfac = 0
  268.           nomid  = 0
  269.           notype = 0
  270.  
  271. C  rho, E, nu pour les massifs
  272.           IF (MFR.EQ.1) THEN
  273.             nbrobl = 3
  274.             SEGINI,nomid
  275.             lesobl(1) = 'RHO '
  276.             lesobl(2) = 'YOUN'
  277.             lesobl(3) = 'NU  '
  278.  
  279.             nbtype = 1
  280.             SEGINI,notype
  281.             type(1) = 'REAL*8'
  282. C
  283. C  rho, cson, rhoref, cref, rlcar pour les liquides
  284.           ELSE IF (MFR.EQ.11.OR.MFR.EQ.41) THEN
  285.             nbrobl = 5
  286.             SEGINI,nomid
  287.             lesobl(1) = 'RHO '
  288.             lesobl(2) = 'CSON'
  289.             lesobl(3) = 'RORF'
  290.             lesobl(4) = 'CREF'
  291.             lesobl(5) = 'LCAR'
  292.  
  293.             nbtype = 1
  294.             SEGINI,notype
  295.             type(1) = 'REAL*8'
  296.           ENDIF
  297.  
  298.           MOMATR = nomid
  299.           MOTYPM = notype
  300.           NMATR = nbrobl
  301.           NMATF = nbrfac
  302.           NMATT = NMATR+NMATF
  303.  
  304. C--- Verification du support des composantes recherchées
  305.           IF (MOMATR.NE.0) THEN
  306.             CALL QUESUQ(IMODEL,IPCHE1,3,0,MOMATR,IPLAZ,ISUPM,iret)
  307.             IF (ISUPM.GT.1) GOTO 1199
  308.           ENDIF
  309. C
  310. C-- Recuperation des noms de composantes CARACTERISTIQUES
  311.           nbrobl = 0
  312.           nbrfac = 0
  313.           nomid  = 0
  314.           notype = 0
  315.  
  316. C Epaisseur du massif en contraintes planes
  317.           IF (MFR.EQ.1 .AND. IFOUR.EQ.-2) THEN
  318.             nbrfac = 1
  319.             SEGINI,nomid
  320.             lesfac(1) = 'DIM3'
  321.             nbtype = 1
  322.             SEGINI,notype
  323.             type(1) = 'REAL*8'
  324.           ENDIF
  325.  
  326.           MOCARA = nomid
  327.           MOTYPC = notype
  328.           NCARA = nbrobl
  329.           NCARF = nbrfac
  330.           NCARR = NCARA+NCARF
  331.  
  332. C--- Verification du support des composantes recherchées
  333.           IF (MOCARA.NE.0) THEN
  334.             CALL QUESUQ(IMODEL,IPCHE1,3,0,MOCARA,IPLAZ,ISUPC,iret)
  335.             IF (ISUPC.GT.1) GOTO 1199
  336.           ENDIF
  337.  
  338. C-- Segment d'integration MINTE
  339.           MINTE = IPMINT
  340.           SEGACT,MINTE
  341.           NBPGAU = POIGAU(/1)
  342.  
  343. C- Partionnement si necessaire de la matrice d'amortissement
  344. C- determinant ainsi le nombre d'objets elementaires de MRIGID
  345.           LTRK = oooval(1,4)
  346.           IF (LTRK.EQ.0) LTRK = oooval(1,1)
  347.       LTRK=MAX(LTRK,2**24)
  348. * Ajout a la taille en mots de la matrice des infos du segment
  349.           LSEG = LRE*LRE*NBELE2 + 16
  350.           NBLPRT = (LSEG-1)/LTRK + 1
  351.           NBLMAX = (NBELE2-1)/NBLPRT + 1
  352.           NBLPRT = (NBELE2-1)/NBLMAX + 1
  353. c*        write(ioimp,*) ' frvisq : nblprt nblmax = ',nblprt,nblmax,nbele2
  354. C*OF : Pour l'instant pas de partition pour FRVISQ
  355.           NBLPRT = 1
  356.  
  357. C-- Ajout de la matrice d'AMORTISSEMENT a la matrice globale
  358.           NRIGE0 = IRIGEL(/2)
  359.           NRIGEL = NRIGE0 + NBLPRT
  360.           SEGADJ,MRIGID
  361.  
  362.           descr  = IDESCR
  363.           meleme = IPOGEO
  364.           nbnn   = NBNOE2
  365.           nbelem = NBELE2
  366.           nbsous = 0
  367.           nbref  = 0
  368.  
  369.           DO 120 irige = 1, NBLPRT
  370.  
  371. C-- Mettre ici la partition du maillage IPOGEO
  372.             ipmail = meleme
  373.             ipdesc = descr
  374.  
  375. C- Initialisation de la matrice de rigidite elementaire (xmatri)
  376.             NELRIG = nbelem
  377.             RIGREL=0
  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.  
  574.  
  575.  

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