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    CB215821  25/09/18    21:15:02     12364          
  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.             SEGINI,xmatri
  378.             ipmatr = xmatri
  379.  
  380. C- Recuperation des valeurs des proprietes materiau et geometriques
  381. c* Note : les proprietes sont les valeurs au support des EF massifs
  382. c*        et non celles au niveau de l'enveloppe surfacique !
  383. c*        Cela ne marche que si les proprietes sont constantes. Dans
  384. c*        les autres cas, le resultat est... Pour eviter cela, on met
  385. c*        un test sur la constance du champ !
  386.             ivamat = 0
  387.             ivacar = 0
  388.             IF (MOMATR.NE.0) THEN
  389.               CALL KOMCHA(IPCHE1,IPT1,CONM,MOMATR,MOTYPM,1,
  390. c*              CALL KOMCHA(IPCHE1,ipmail,CONM,MOMATR,MOTYPM,1,
  391.      &                    INFOS,3, ivamat)
  392.               IF (IERR.NE.0) GOTO 1199
  393.               mptval = ivamat
  394.               do i = 1, NMATT
  395.                 if (ival(i).ne.0) then
  396.                   melval = IVAL(i)
  397.                   if (velche(/1).ne.1 .and. velche(/2).ne.1) then
  398.                     write(ioimp,*) 'Champ MATERIAU non constant'
  399.                     call erreur(21)
  400.                     goto 1199
  401.                   endif
  402.                 endif
  403.               enddo
  404.               IF (ISUPM.EQ.1) THEN
  405.                 CALL VALCHE(ivamat,NMATT,IPMINT,IPPORE,MOMATR,MELE)
  406.                 IF (IERR.NE.0) THEN
  407.                   ISUPM = 0
  408.                   GOTO 1199
  409.                 ENDIF
  410.               ENDIF
  411.             ENDIF
  412.             IF (MOCARA.NE.0) THEN
  413.               CALL KOMCHA(IPCHE1,IPT1,CONM,MOCARA,MOTYPC,1,
  414. c*              CALL KOMCHA(IPCHE1,ipmail,CONM,MOCARA,MOTYPC,1,
  415.      &                    INFOS,3, ivacar)
  416.               IF (IERR.NE.0) GOTO 1199
  417.               mptval = ivacar
  418.               do i = 1, NCARR
  419.                 if (ival(i).ne.0) then
  420.                   melval = IVAL(i)
  421.                   if (velche(/1).ne.1 .and. velche(/2).ne.1) then
  422.                     write(ioimp,*) 'Champ MATERIAU non constant'
  423.                     call erreur(21)
  424.                     goto 1199
  425.                   endif
  426.                 endif
  427.               enddo
  428.               IF (ISUPC.EQ.1) THEN
  429.                 CALL VALCHE(ivacar,NCARR,IPMINT,IPPORE,MOCARA,MELE)
  430.                 IF (IERR.NE.0)THEN
  431.                   ISUPC = 0
  432.                   GOTO 1199
  433.                 ENDIF
  434.               ENDIF
  435.             ENDIF
  436.  
  437. C  distinction des cas 2D et 3D
  438. C______________________________________________________________________C
  439. C                                                                      C
  440. C            CAS DES ELEMENTS MASSIFS BIDIMENSIONNELS                  C
  441. C                  FACES ASSOCIEES SEG2 OU SEG3                        C
  442. C______________________________________________________________________C
  443. C                                                                      C
  444.          IF (MELE.EQ.2.OR.MELE.EQ.3) THEN
  445. C
  446.             CALL FROA2D(ipmail,ipmatr,IPMINT,ivamat,ivacar,
  447.      1                  MELE,MFR,LRE,NDDL)
  448. C
  449. C______________________________________________________________________C
  450. C                                                                      C
  451. C               CAS DES ELEMENTS LIQUIDES 2D OU 3D                     C
  452. C               FACES ASSOCIEES LSE2, LTR3 OU LQU4                     C
  453. C______________________________________________________________________C
  454. C                                                                      C
  455. C
  456.          ELSE IF(MELE.EQ.97.OR.MELE.EQ.35.OR.MELE.EQ.36) THEN
  457. C
  458.             CALL LFROA(ipmail,ipmatr,IPMINT,ivamat,ivacar,
  459.      1                  MELE,MFR,LRE,NDDL)
  460. C
  461. C______________________________________________________________________C
  462. C                                                                      C
  463. C            CAS DES ELEMENTS MASSIFS TRIDIMENSIONNELS                 C
  464. C             FACES ASSOCIEES FAC3,FAC4,FAC6 OU FAC8                   C
  465. C______________________________________________________________________C
  466. C
  467.          ELSE IF(MELE.EQ.31.OR.MELE.EQ.32.OR.MELE.EQ.33.
  468.      1                                   OR.MELE.EQ.34)THEN
  469. C
  470.             CALL FROA3D(ipmail,ipmatr,IPMINT,ivamat,ivacar,
  471.      1                  MELE,MFR,LRE,NDDL)
  472. C
  473. C  erreur, l'élément n'est pas encore implémenté
  474. C
  475.          ELSE
  476. C
  477.             MOTERR(1:4)=NOMTP(MELE)
  478.             MOTERR(5:12)='FRVISQ'
  479.             CALL ERREUR (86)
  480.          ENDIF
  481. C
  482.          IF (ISUPM.EQ.1 .OR. nblprt.GT.1) THEN
  483.            CALL DTMVAL(ivamat,3)
  484.          ELSE
  485.            CALL DTMVAL(ivamat,1)
  486.          ENDIF
  487.          IF (ISUPC.EQ.1 .OR. nblprt.GT.1) THEN
  488.            CALL DTMVAL(ivacar,3)
  489.          ELSE
  490.            CALL DTMVAL(ivacar,1)
  491.          ENDIF
  492.  
  493.           xmatri = ipmatr
  494.           IF (NBLPRT.GT.1) THEN
  495.             meleme = ipmail
  496.           ENDIF
  497.  
  498. C- Sortie prematuree en cas d'erreur
  499.           IF (IERR.NE.0) GOTO 1199
  500.  
  501. C- Stockage de la matrice
  502.           jrige = NRIGE0 + irige
  503.           COERIG(jrige)   = 1.
  504.           IRIGEL(1,jrige) = ipmail
  505.           IRIGEL(2,jrige) = 0
  506.           IRIGEL(3,jrige) = ipdesc
  507.           IRIGEL(4,jrige) = ipmatr
  508.           IRIGEL(5,jrige) = NIFOUR
  509.           IRIGEL(6,jrige) = 0
  510.           IRIGEL(7,jrige) = 0
  511. * matrice non symetrique (forces sur pi seulement
  512. * qui dependent de p)
  513.           IF (MFR.EQ.11.OR.MFR.EQ.41) THEN
  514.             IRIGEL(7,jrige) = 2
  515.             xmatri.symre=2
  516.            ENDIF
  517.           SEGDES,xmatri
  518.           IRIGEL(8,jrige) = 0
  519.  
  520.  120      CONTINUE
  521. C- Fin de la boucle de partition maillage/rigidite
  522.  
  523.  1199     CONTINUE
  524.           IF (MOMATR.NE.0) THEN
  525.             nomid = MOMATR
  526.             SEGSUP,nomid
  527.             notype = MOTYPM
  528.             SEGSUP,notype
  529.           ENDIF
  530.           IF (MOCARA.NE.0) THEN
  531.             nomid = MOCARA
  532.             SEGSUP,nomid
  533.             notype = MOTYPC
  534.             SEGSUP,notype
  535.           ENDIF
  536. C
  537.  
  538. C- Sortie prematuree en cas d'erreur
  539.           IF (IERR.NE.0) GOTO 1098
  540.  
  541.  110    CONTINUE
  542. C- Fin de la boucle sur (les sous-zones de) l'enveloppe
  543. C
  544.  1098   CONTINUE
  545.  1099   CONTINUE
  546. C- Sortie prematuree en cas d'erreur
  547.         IF (IERR.NE.0) GOTO 999
  548. C
  549.  100  CONTINUE
  550.       SEGSUP,IMODE1
  551. C- Fin de la boucle sur les modeles elementaires
  552. C
  553.       NRIGE0 = IRIGEL(/2)
  554.       IF (NRIGE0.EQ.0) THEN
  555.         CALL ERREUR(902)
  556.       ENDIF
  557.  
  558.  999  CONTINUE
  559.       IF (IERR.EQ.0) THEN
  560.         IPRIG = MRIGID
  561.         SEGDES,MRIGID
  562.       ELSE
  563.         IPRIG = 0
  564.         SEGSUP,MRIGID
  565.       ENDIF
  566.  
  567. c      RETURN
  568.       END
  569.  
  570.  
  571.  
  572.  
  573.  

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