C FPMA3D    SOURCE    OF166741  25/02/21    21:16:43     12166          

C____________________________________________________________________
C   CALCULE LES FORCES DE PRESSIONS SUR LES FACES D ELEMENTS
C   MASSIFS TRIDIMENSIONNELS
C
C   ENTREES :
C   ---------
C
C    IPTVPR  POINTEUR SUR UN MELVAL CONTENANT LES PRESSIONS APPLIQUEES
C            0 SI ON A DONNE UNE PRESSION CONSTANTE
C    IPMAIL  POINTEUR SUR UN OBJET GEOMETRIQUE
C    IPTINT  POINTEUR SUR UN MINTE CONTENANT LES POINTS D INTEGRATION
C            ACTIF EN ENTREE ET EN SORTIE SANS MODIFICATION
C    IVAFOR  POINTEUR SUR UN MPTVAL ET LES MELVAL CONTENANT LES FORCES
C            NODALES RESUL
C
C      JACQUELINE BROCHARD  AVRIL 85
C
C      PASSAGE AUX NOUVEAUX CHAMELEM PAR JM CAMPENON LE 17 09 90
C
C______________________________________________________________________

      SUBROUTINE FPMA3D(IPTVPR,IPMAIL,ipmaim,IPTINT,IVAFOR,XP
     &                 ,netn1,ietn1)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

-INC PPARAM
-INC CCOPTIO

-INC SMCHAML
-INC SMELEME
-INC SMINTE
-INC SMCOORD

-INC TMPTVAL

      segment netn(notn)
      segment ietn(letn)

      SEGMENT WORK
        REAL*8 XE(3,NBNN)
      ENDSEGMENT

      idimp1 = IDIM+1
*   prob optimiseur il faut initialiser melva1
      MELVA1=IVAFOR
      IF (IPTVPR.NE.0) THEN
        MELVA1=IPTVPR
c*        SEGACT,MELVA1                            <- ACTIF EN E/S
        IG11 = MELVA1.VELCHE(/1)
        IB12 = MELVA1.VELCHE(/2)
      ENDIF

      MINTE=IPTINT
C*      SEGACT,MINTE                               <- ACTIF EN E/S
      NBPGAU=POIGAU(/1)

      MELEME = IPMAIL
c*      SEGACT,MELEME                              <- ACTIF EN E/S
      NBNN   = meleme.NUM(/1)
      NBELEM = meleme.NUM(/2)

      SEGINI,WORK

      netn = netn1
      ietn = ietn1
      ipt1 = ipmaim
      IF (IPT1.GT.0) THEN
        if (netn.eq.0 .or. ietn.eq.0) then
          write(ioimp,*) 'FPMA3D : incompatibilite netn, ietn & IPMAIM'
        endif
c*        SEGACT,IPT1                                <- ACTIF en E/S
        nbnn1 = ipt1.num(/1)
        nbel1 = ipt1.num(/2)
      ELSE
        if (netn.gt.0 .or. ietn.gt.0) then
          write(ioimp,*) 'FPMA2D : incompatibilite netn, ietn & IPMAIM'
        endif
      ENDIF
C
C     BOUCLE SUR LES ELEMENTS
C
      DO IB = 1, NBELEM

        CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)

        XFLOT = +1.D0
        IF (netn.GT.0) THEN
          DO inf = 1, nbnn
            ip = meleme.num(inf,ib)
            ideb = netn(ip)+1
            ifin = netn(ip+1)
            do itn = ideb, ifin
              IEM = ietn(itn)
              jne = 0
              do i = 1, nbnn
                ino = num(i,ib)
                do i1 = 1, nbnn1
                  if (ino.eq.ipt1.num(i1,IEM)) jne=jne+1
                enddo
              enddo
              if (jne.eq.nbnn) goto 170
            enddo
          ENDDO
          CALL ERREUR(26)
          GOTO 9900
 170      continue
          XG = 0.D0
          YG = 0.D0
          ZG = 0.D0
          DO I = 1, NBNN1
            ino = (IPT1.NUM(I,IEM)-1)*idimp1
            XG=XG+XCOOR(ino+1)
            YG=YG+XCOOR(ino+2)
            ZG=ZG+XCOOR(ino+3)
          ENDDO
          XG=XG / NBNN1
          YG=YG / NBNN1
          ZG=ZG / NBNN1

          XK=0.D0
          YK=0.D0
          ZK=0.D0
          DO i = 1,NBNN
            XK=XK+XE(1,I)
            YK=YK+XE(2,I)
            ZK=ZK+XE(3,I)
          ENDDO
          XK=XK/NBNN
          YK=YK/NBNN
          ZK=ZK/NBNN

          V_1 = XG - XK
          V_2 = YG - YK
          V_3 = ZG - ZK
          r_z = 1.D0 / SQRT(V_1*V_1+V_2*V_2+V_3*V_3)
          V_1 = V_1 * r_z
          V_2 = V_2 * r_z
          V_3 = V_3 * r_z
        ENDIF
C
C     BOUCLE SUR LES POINTS DE GAUSS
C
        DO IGAU = 1, NBPGAU

          VNQSI1 = 0.D0
          VNQSI2 = 0.D0
          VNQSI3 = 0.D0
          VNETA1 = 0.D0
          VNETA2 = 0.D0
          VNETA3 = 0.D0
          DO I = 1, NBNN
            VNQSI1 = VNQSI1+SHPTOT(2,I,IGAU)*XE(1,I)
            VNQSI2 = VNQSI2+SHPTOT(2,I,IGAU)*XE(2,I)
            VNQSI3 = VNQSI3+SHPTOT(2,I,IGAU)*XE(3,I)
            VNETA1 = VNETA1+SHPTOT(3,I,IGAU)*XE(1,I)
            VNETA2 = VNETA2+SHPTOT(3,I,IGAU)*XE(2,I)
            VNETA3 = VNETA3+SHPTOT(3,I,IGAU)*XE(3,I)
          ENDDO

          VNN1 = VNQSI2*VNETA3-VNQSI3*VNETA2
          VNN2 = VNQSI3*VNETA1-VNQSI1*VNETA3
          VNN3 = VNQSI1*VNETA2-VNQSI2*VNETA1

          if (igau.eq.1.and.netn.gt.0) then
            vnnn = 1.D0 / SQRT(vnn1*vnn1+vnn2*vnn2+vnn3*vnn3)
            test = v_1*(vnn1*vnnn) + v_2*(vnn2*vnnn) + v_3*(vnn3*vnnn)
            if (test.lt.0d0) xflot = -1.d0
          endif

          r_z = POIGAU(IGAU) * XFLOT
          IF (IPTVPR.NE.0) THEN
            IGMN=MIN(IGAU,IG11)
            IBMN=MIN(IB  ,IB12)
            r_z = r_z * MELVA1.VELCHE(IGMN,IBMN)
          ELSE
            r_z = r_z * XP
          ENDIF

          T1 = r_z * VNN1
          T2 = r_z * VNN2
          T3 = r_z * VNN3

          MPTVAL=IVAFOR
          MELVAL=IVAL(1)
          DO i=1,NBNN
            VELCHE(i,IB) = VELCHE(i,IB) + SHPTOT(1,i,IGAU)*T1
          ENDDO
          MELVAL=IVAL(2)
          DO i=1,NBNN
            VELCHE(i,IB) = VELCHE(i,IB) + SHPTOT(1,i,IGAU)*T2
          ENDDO
          MELVAL=IVAL(3)
          DO i=1,NBNN
            VELCHE(i,IB) = VELCHE(i,IB) + SHPTOT(1,i,IGAU)*T3
          ENDDO

        ENDDO

      ENDDO

 9900 CONTINUE
      SEGSUP,WORK

c*      SEGDES,MINTE                               <- ACTIF en E/S
c*      SEGDES,MELEME                              <- ACTIF en E/S
c*      IF (IPTVPR.NE.0) SEGDES,MELVA1             <- ACTIF en E/S

      RETURN
      END