Sources Esope | Cast3M

lfroa
C LFROA     SOURCE    OF166741  25/02/21    21:17:52     12166          
 
      SUBROUTINE LFROA(IPOGEO,IPMATR,IPMINT,IVAMAT,
     1                  IVACAR,MELE,MFR,LRE,NDDL)
C
C***********************************************************************
C                                                                      *
C  Routine appelée par FRVISQ.                                         *
C                                                                      *
C  Calcule l'amortissement de frontière pour les liquides de face LSE2 *
C  dans le cas 2D, ou de face LTR3 ou LQU4 dans le cas 3D.             *
C                                                                      *
C  Entrées :                                                           *
C  --------                                                            *
C                                                                      *
C     IPOGEO : pointeur sur le maillage de l'enveloppe des massifs,    *
C              type MELEME                                             *
C     IPMATR : pointeur sur le segment IMATRI, chapeau des rigidités   *
C              élémentaires                                            *
C     IPMINT : pointeur sur le segment d'intégration, type MINTE       *
C     IVAMAT : pointeur sur un segment MPTVAL de données matériau      *
C     IVACAR : pointeur sur un segment MPTVAL de caractéristiques      *
C              (épaisseur dans le cas contraintes planes)              *
C     MELE   : numéro de l'élément fini associé à la face du massif    *
C     MFR    : numéro de la formulation                                *
C     LRE    : taille de la matrice d'amortissement à construire       *
C     NDDL   : nombre de degrés de liberté                             *
C                                                                      *
C  Remplit le segment XMATRI pour chaque élément de la sous-zone.      *
C***********************************************************************
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
 
-INC SMRIGID
-INC SMELEME
-INC SMCOORD
-INC SMCHAML
-INC SMINTE
 
-INC TMPTVAL
 
      SEGMENT,MWORK
        REAL*8 XE(3,NBNN)
        REAL*8 REL(LRE,LRE)
        REAL*8 SHPWRK(6,NBNN),BGENE(NDDL,LRE)
        REAL*8 VALMAT(NV1)
        REAL*8 RCLMAT(NDDL,NDDL)
      ENDSEGMENT
 
      MELEME=IPOGEO
c*      SEGACT MELEME
      NBNN=NUM(/1)
      NBELEM=NUM(/2)
C
      MINTE=IPMINT
c*      SEGACT,MINTE
      NBPGAU=POIGAU(/1)
C
      DIM3=1.D0
C
      xMATRI=IPMATR
c*      SEGACT,xMATRI*MOD
c*      NLIGRD=LRE
c*      NLIGRP=LRE
 
      NV1=5
      SEGINI,MWORK
C
C  boucle sur les éléments
C
      DO 1 IB=1,NBELEM
C
C  on cherche les coordonnées de l'élément IB
C
        CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
        CALL ZERO(REL,LRE,LRE)
        CALL ZERO(RCLMAT,NDDL,NDDL)
C
C  boucle sur les points de Gauss
C
        DO 10 IGAU=1,NBPGAU
C
C  récupération de l'épaisseur
C
          IF (IFOUR.EQ.-2) THEN
            IF (IVACAR.NE.0) THEN
              MPTVAL=IVACAR
              MELVAL=IVAL(1)
              IF (MELVAL.NE.0) THEN
                IGMN=MIN(IGAU,VELCHE(/1))
                IBMN=MIN(IB,VELCHE(/2))
                DIM3=VELCHE(IGMN,IBMN)
              ELSE
                DIM3=1.D0
              ENDIF
            ENDIF
          ENDIF
C
C  récupération des données matériau
C
         MPTVAL=IVAMAT
         DO 11 J=1,5
            MELVAL=IVAL(J)
            IGMN=MIN(IGAU,VELCHE(/1))
            IBMN=MIN(IB,VELCHE(/2))
            VALMAT(J)=VELCHE(IGMN,IBMN)
 11      CONTINUE
C
         RHO=VALMAT(1)
         CSON=VALMAT(2)
         ROREF=VALMAT(3)
         CREF=VALMAT(4)
         RLCAR=VALMAT(5)
C
C  coefficient d'amortissement normalisé
C
         RCL=(1.D0/CSON/rho)*ROREF*RLCAR*roref*cref**2/rlcar
C
C  matrice RCLMAT
C
         RCLMAT(1,1)=RCL
         RCLMAT(2,2)=RCL
C
C  calcul de la matrice N des fonctions de forme
C
         XDPGE=0.D0
         YDPGE=0.D0
         CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL,
     1   DIM3,XE,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE)
C
C  mise à zéro des composantes sur p
C
            DO 21 J=1,LRE
               BGENE(2,J)=0.D0
 21         CONTINUE
C
C  calcul du jacobien
C
          IF (IDIM.EQ.3) THEN
            DXDQSI=0.D0
            DYDQSI=0.D0
            DZDQSI=0.D0
            DXDETA=0.D0
            DYDETA=0.D0
            DZDETA=0.D0
            DO 40 I=1,NBNN
              DXDQSI=DXDQSI+SHPTOT(2,I,IGAU)*XE(1,I)
              DYDQSI=DYDQSI+SHPTOT(2,I,IGAU)*XE(2,I)
              DXDETA=DXDETA+SHPTOT(3,I,IGAU)*XE(1,I)
              DYDETA=DYDETA+SHPTOT(3,I,IGAU)*XE(2,I)
              DZDQSI=DZDQSI+SHPTOT(2,I,IGAU)*XE(3,I)
              DZDETA=DZDETA+SHPTOT(3,I,IGAU)*XE(3,I)
 40         CONTINUE
            DJAC=SQRT((DYDQSI*DZDETA-DYDETA*DZDQSI)**2+
     1                (DXDETA*DZDQSI-DXDQSI*DZDETA)**2+
     2                (DXDQSI*DYDETA-DXDETA*DYDQSI)**2)
          ELSE
c*          ELSE IF(IDIM.EQ.2) THEN
            DXDQSI=0.D0
            DYDQSI=0.D0
            DO I=1,NBNN
              DXDQSI=DXDQSI+SHPTOT(2,I,IGAU)*XE(1,I)
              DYDQSI=DYDQSI+SHPTOT(2,I,IGAU)*XE(2,I)
            ENDDO
            DJAC=SQRT(DXDQSI**2+DYDQSI**2)
          ENDIF
C
C  calcul de l'élément de volume
C
         IF (IFOUR.LT.0.OR.IFOUR.EQ.2) THEN
           R=1.D0
           IF (IFOUR.EQ.-2) R=DIM3
         ELSE
           R=0.D0
           DO I=1,NBNN
              R=R+SHPTOT(1,I,IGAU)*XE(1,I)
           ENDDO
           IF (IFOUR.EQ.0.OR.(IFOUR.EQ.1
     1              .AND.NIFOUR.EQ.0)) THEN
               R=2*XPI*R
           ELSE
c*           ELSEIF (IFOUR.EQ.1.AND.NIFOUR.NE.0) THEN
               R=XPI*R
           ENDIF
         ENDIF
C
C  construction de la matrice d'amortissement
C
         DJAC=DJAC*POIGAU(IGAU)*R
         CALL BDBST(BGENE,DJAC,RCLMAT,LRE,NDDL,REL)
C
 10   CONTINUE
C
C  remplissage de XMATRI
C
      CALL REMPMT(REL,LRE,RE(1,1,ib))
*
** la matrice calculée ci dessus serait correcte si l'on avait
** une formulation uniquement en pi. Comme on retient la formulation
** en p et pi et pour ne rien ajouter sur la ligne en pi (si non on
** ne satisfait plus la relation entre p et pi) on aboutit à une matrice
** disymetrique avec des termes uniquement sur le lignes en pi et collones
** en p. D'où la matrice suivante :
 
       do 30 i = 1, lre
       do 31 j = 2, lre
         ix = mod(j,2)
        if (ix.eq.0) then
         re(i,j,ib) = -re( i,(j-1),ib)
         re( i,(j-1),ib) = 0.d0
        endif
 31   continue
 30   continue
 
 1    CONTINUE
 
      SEGSUP,MWORK
c*      SEGDES MELEME,MINTE,xMATRI
 
      RETURN
      END