mass-fluid_qd.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    mass_fluid_qd   ../starter/source/fluid/mass-fluid_qd.F
!||--- called by ------------------------------------------------------
!||    hm_read_bem     ../starter/source/loads/bem/hm_read_bem.F
!||--- calls      -----------------------------------------------------
!||    ancmsg          ../starter/source/output/message/message.F
!||    intanl_qd       ../starter/source/fluid/intanl_qd.F
!||    intanl_tg       ../starter/source/fluid/intanl_tg.F
!||    intgqd          ../starter/source/fluid/inte_qd.F
!||    intgtg          ../starter/source/fluid/inte_tg.F
!||    inthqd          ../starter/source/fluid/inte_qd.F
!||    inthtg          ../starter/source/fluid/inte_tg.F
!||    invert          ../starter/source/constraints/general/rbe3/hm_read_rbe3.F
!||--- uses       -----------------------------------------------------
!||    message_mod     ../starter/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE mass_fluid_qd(NNO, NEL, IFLOW, IBUF, ELEM, X,
     .                         NORMAL, AF, MFLE, CBEM, RHO,IRESP)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "units_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NNO, NEL, IFLOW(*), IBUF(*), ELEM(5,*)
      my_real x(3,*), af(*), normal(3,*), mfle(nel,*), cbem(nel,*), rho
      INTEGER,INTENT(IN):: IRESP
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER IFORM, KFORM, ILVOUT, NELMAX
      INTEGER IN, JN, KN, N1, N2, N3, N4, N5, NN1, NN2, NN3, NN4, NNJ, IEL, JEL, ERR
      my_real x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, xq, yq, zq, r2,
     .        xp, yp, zp, x13, y13, z13, x24, y24, z24,
     .        nrx, nry, nrz,  area, d2, rval, rvlh, rvlg
      my_real massx, massy, massz, wi(4,2), sum
      my_real, DIMENSION(:,:), ALLOCATABLE :: bbem, ebem
      my_real :: alpha, beta
      INTEGER :: AERR
      iform  = iflow(13)
      ilvout = iflow(17)
      kform  = iflow(23)
      nelmax = iflow(28)
C-------------------------------------------------------------------
C 1- Compute Bij Cij and Ai
C-------------------------------------------------------------------
      IF (ilvout>=1) THEN
       WRITE(istdo,'(A)') ' .. FLUID MASS MATRIX:
     . ASSEMBLY OF INTEGRAL OPERATORS'
      ENDIF
C Collocation BEM 
      wi(1,1)=fourth
      wi(2,1)=fourth
      wi(3,1)=fourth
      wi(4,1)=fourth
      wi(1,2)=third
      wi(2,2)=third
      wi(3,2)=one_over_6
      wi(4,2)=one_over_6
      ALLOCATE(bbem(nel,nel), stat = aerr)
      IF (aerr /= 0) THEN
         CALL ancmsg(msgid = 1710, anmode=aninfo, msgtype = msgerror) 
      ENDIF
      DO in=1,nel
         DO jn=1,nel
            bbem(in,jn)=zero
            cbem(in,jn)=zero
         ENDDO
      ENDDO
 
      IF(iform == 1) THEN
C Gauss integration
        DO iel=1,nel
C           IF (ILVOUT>=2) CALL PROGBAR_C(IEL,NEL)
           n1=elem(1,iel)
           n2=elem(2,iel)
           n3=elem(3,iel)
           n4=elem(4,iel)
           n5=elem(5,iel)
           nn1=ibuf(n1)
           nn2=ibuf(n2)
           nn3=ibuf(n3)
           nn4=ibuf(n4)
           x1=x(1,nn1)
           x2=x(1,nn2)
           x3=x(1,nn3)
           x4=x(1,nn4)
           y1=x(2,nn1)
           y2=x(2,nn2)
           y3=x(2,nn3)
           y4=x(2,nn4)
           z1=x(3,nn1)
           z2=x(3,nn2)
           z3=x(3,nn3)
           z4=x(3,nn4)
C  Control point P
           xp=wi(1,n5)*x1+wi(2,n5)*x2+wi(3,n5)*x3+wi(4,n5)*x4
           yp=wi(1,n5)*y1+wi(2,n5)*y2+wi(3,n5)*y3+wi(4,n5)*y4
           zp=wi(1,n5)*z1+wi(2,n5)*z2+wi(3,n5)*z3+wi(4,n5)*z4
           d2=min((xp-x1)**2+(yp-y1)**2+(zp-z1)**2,
     .            (xp-x2)**2+(yp-y2)**2+(zp-z2)**2,
     .            (xp-x3)**2+(yp-y3)**2+(zp-z3)**2,
     .            (xp-x4)**2+(yp-y4)**2+(zp-z4)**2)
           nrx=normal(1,iel)
           nry=normal(2,iel)
           nrz=normal(3,iel)
C   
           DO jel=1,nel
              IF(iel == jel) THEN
                bbem(iel,jel)=two*sqrt(pi*af(jel))
                cbem(iel,jel)=two*pi
              ELSE
                n1=elem(1,jel)
                n2=elem(2,jel)
                n3=elem(3,jel)
                n4=elem(4,jel)
                jn=elem(5,jel)
                nn1=ibuf(n1)
                nn2=ibuf(n2)
                nn3=ibuf(n3)
                nn4=ibuf(n4)
                x1=x(1,nn1)
                x2=x(1,nn2)
                x3=x(1,nn3)
                x4=x(1,nn4)
                y1=x(2,nn1)
                y2=x(2,nn2)
                y3=x(2,nn3)
                y4=x(2,nn4)
                z1=x(3,nn1)
                z2=x(3,nn2)
                z3=x(3,nn3)
                z4=x(3,nn4)
                area=af(jel)
C  Position de la source Q
                xq=wi(1,jn)*x1+wi(2,jn)*x2+wi(3,jn)*x3+wi(4,jn)*x4
                yq=wi(1,jn)*y1+wi(2,jn)*y2+wi(3,jn)*y3+wi(4,jn)*y4
                zq=wi(1,jn)*z1+wi(2,jn)*z2+wi(3,jn)*z3+wi(4,jn)*z4
 
                IF(jn == 1) THEN
                  CALL inthqd(x1 , y1,  z1,  x2,  y2,  z2,
     .                        x3,  y3,  z3,  x4,  y4,  z4,
     .                        xp,  yp,  zp,  xq,  yq,  zq,
     .                        nrx, nry, nrz, d2, area, rval)
                  cbem(iel,jel)=rval
                  CALL intgqd(x1 , y1,  z1,  x2,  y2,  z2,
     .                        x3,  y3,  z3,  x4,  y4,  z4,
     .                        xp,  yp,  zp,  xq,  yq,  zq,
     .                        d2, area, rval)
                  bbem(iel,jel)=rval
                ELSEIF(jn == 2) THEN
                  CALL inthtg(x1 , y1,  z1,  x2,  y2,  z2,  
     .                        x3,  y3,  z3,  xp,  yp,  zp,                    
     .                        nrx,nry, nrz,  d2,  area,
     .                        xq,  yq,  zq, rval )
                  cbem(iel,jel)=rval
                  CALL intgtg(x1,  y1,  z1,  x2,  y2,  z2,  
     .                        x3,  y3,  z3,  xp,  yp,  zp,                    
     .                        d2,  area,
     .                        xq,  yq,  zq, rval )
                  bbem(iel,jel)=rval
                ENDIF
              ENDIF
           ENDDO
        ENDDO 
 
      ELSEIF(iform == 2) THEN
C Integration analytique
        DO iel=1,nel
C           IF (ILVOUT>=2) CALL PROGBAR_C(IEL,NEL)
           n1=elem(1,iel)
           n2=elem(2,iel)
           n3=elem(3,iel)
           n4=elem(4,iel)
           in=elem(5,iel)
           nn1=ibuf(n1)
           nn2=ibuf(n2)
           nn3=ibuf(n3)
           nn4=ibuf(n4)
           x1=x(1,nn1)
           x2=x(1,nn2)
           x3=x(1,nn3)
           x4=x(1,nn4)
           y1=x(2,nn1)
           y2=x(2,nn2)
           y3=x(2,nn3)
           y4=x(2,nn4)
           z1=x(3,nn1)
           z2=x(3,nn2)
           z3=x(3,nn3)
           z4=x(3,nn4)
C  Control point P
           xp=wi(1,in)*x1+wi(2,in)*x2+wi(3,in)*x3+wi(4,in)*x4
           yp=wi(1,in)*y1+wi(2,in)*y2+wi(3,in)*y3+wi(4,in)*y4
           zp=wi(1,in)*z1+wi(2,in)*z2+wi(3,in)*z3+wi(4,in)*z4
           DO jel=1,nel
              n1=elem(1,jel)
              n2=elem(2,jel)
              n3=elem(3,jel)
              n4=elem(4,jel)
              jn=elem(5,jel)
              nn1=ibuf(n1)
              nn2=ibuf(n2)
              nn3=ibuf(n3)
              nn4=ibuf(n4)
              x1=x(1,nn1)
              x2=x(1,nn2)
              x3=x(1,nn3)
              x4=x(1,nn4)
              y1=x(2,nn1)
              y2=x(2,nn2)
              y3=x(2,nn3)
              y4=x(2,nn4)
              z1=x(3,nn1)
              z2=x(3,nn2)
              z3=x(3,nn3)
              z4=x(3,nn4)
              nrx=normal(1,jel)
              nry=normal(2,jel)
              nrz=normal(3,jel)
              area=af(jel)
              IF(jn == 1) THEN
                CALL intanl_qd(x1, y1, z1, x2, y2, z2, x3, y3, z3,
     .                         x4, y4, z4, xp, yp, zp, nrx,nry,nrz,
     .                         area, rvlh, rvlg, iel, jel )
 
              ELSEIF(jn == 2) THEN
                CALL intanl_tg(x1, y1, z1, x2, y2, z2, x3, y3, z3,
     .                         xp, yp, zp, nrx,nry,nrz,                     
     .                         area, rvlh, rvlg, iel, jel )
              ENDIF
C Matrice C=H
              IF(iel == jel) THEN
                 cbem(iel,jel)=two*pi
              ELSE
                 cbem(iel,jel)=rvlh
              ENDIF
C Matrice B=G
              bbem(iel,jel)=rvlg
           ENDDO
        ENDDO
      ENDIF
 
      IF(ilvout>=3) THEN
        WRITE (*,'(//,A)') ' BBEM MATRIX'
        IF(nel < 11) THEN
          DO in = 1,nel
             WRITE (*,'(10E13.5)')  (bbem(in,jn),jn=1,nel)
          ENDDO
        ELSE
          DO in = 1,10
             WRITE (*,'(10E13.5)')  (bbem(in,jn),jn=1,10)
          ENDDO
          WRITE (*,'(//,A)') ' BBEM MATRIX B1I'
          WRITE (*,'(10E13.5)')  (bbem(1,jn),jn=1,nel)
        ENDIF
 
        WRITE (*,'(//,A)') ' CBEM MATRIX'
        IF(nel < 11) THEN
          DO in = 1,nel
             WRITE (*,'(10E13.5)')  (cbem(in,jn),jn=1,nel)
          ENDDO
        ELSE
          DO in = 1,10
             WRITE (*,'(10e13.5)')  (CBEM(IN,JN),JN=1,10)
          ENDDO
          WRITE (*,'(//,a)') ' cbem matrix c1i'
          WRITE (*,'(10e13.5)')  (CBEM(1,JN),JN=1,NEL)
       ENDIF
      ENDIF
C------------------------------------------------------------------------
C 2- Compute matrix MFLE = CBEMt * BBEM + BBEMt * CBEM (Case DAA KFORM=1)
C------------------------------------------------------------------------
.AND.      IF(KFORM == 1  NEL >= NELMAX) THEN
c         MFLE  = CBEMt*BBEM ! DGEMM C := alpha*op( A )*op( B ) +beta*C
C                              C = MFLE ; beta = 0
C                              op = T A = CBEM
C                              op = N , B = BBEM
c         MFLE = MFLE + BBEMt*CBEM 
C                              C = MFLE ; beta = 1
C                              op = N A = BBEMt
C                              op = T , B = CBEM
 
C         CMAT(1:NEL,1:NEL) = TRANSPOSE(BBEM(1:NEL,1:NEL))
C         MFLE(1:NEL,1:NEL) = MATMUL(CMAT(1:NEL,1:NEL),CBEM(1:NEL,1:NEL))
C         CMAT(1:NEL,1:NEL) = TRANSPOSE(CBEM(1:NEL,1:NEL))
C         MFLE(1:NEL,1:NEL) = MFLE(1:NEL,1:NEL) + MATMUL(CMAT(1:NEL,1:NEL),BBEM(1:NEL,1:NEL))
C         CMAT(1:NEL,1:NEL) = CBEM(1:NEL,1:NEL)
  
         IF (IRESP == 0)THEN
         ! Double Precision version : use either Lapack / MKP or ARMPL
 
           ALPHA = 1.0D0
           BETA = 0.0D0
           CALL DGEMM('t','n',NEL,NEL,NEL,ALPHA,CBEM,NEL,BBEM,NEL,BETA,MFLE,NEL)
           ALPHA = 1.0D0
           BETA = 1.0D0
           CALL DGEMM('n','t',NEL,NEL,NEL,ALPHA,BBEM,NEL,CBEM,NEL,BETA,MFLE,NEL)
         ELSE
         ! Single Precision version / Bad performance 
           DO JN=1,NEL
             DO IN=1,NEL
               SUM=ZERO
               DO KN=1,NEL
                  SUM=SUM+BBEM(KN,IN)*CBEM(KN,JN)+CBEM(KN,IN)*BBEM(KN,JN)
               ENDDO
               MFLE(IN,JN)=SUM
             ENDDO
           ENDDO
         ENDIF
         RETURN
      ENDIF
C---------------------------------------------------
C 2- Compute fluid mass matrix MFLE
C---------------------------------------------------
      IF (ILVOUT>=1) WRITE(ISTDO,'(a)') ' .. fluid mass matrix'
      ALLOCATE(EBEM(NEL,NEL), STAT = AERR)
      IF (AERR /= 0) THEN
         CALL ANCMSG(MSGID = 1710, ANMODE=ANINFO, MSGTYPE = MSGERROR) 
      ENDIF
       CALL INVERT(CBEM,EBEM,NEL,ERR)
 
      IF(ILVOUT>=3) THEN
        WRITE (*,'(//,a)') ' cbem-1 matrix'
        IF(NEL < 11) THEN
          DO IN = 1,NEL
             WRITE (*,'(10e13.5)')  (EBEM(IN,JN),JN=1,NEL)
          END DO
        ELSE
          DO IN = 1,10
             WRITE (*,'(10e13.5)')  (EBEM(IN,JN),JN=1,10)
          END DO
          WRITE (*,'(//,a)') ' cbem-1 matrix c1i'
          WRITE (*,'(10e13.5)')  (EBEM(1,JN),JN=1,NEL)
        ENDIF
      ENDIF
 
      DO IN=1,NEL
         DO JN=1,NEL
            MFLE(IN,JN)=ZERO
            CBEM(IN,JN)=ZERO
         ENDDO
      ENDDO
      CBEM(1:NEL, 1:NEL) = MATMUL(BBEM(1:NEL, 1:NEL), EBEM(1:NEL, 1:NEL))
c$$$      DO IN=1,NEL
c$$$         DO JN=1,NEL
c$$$            DO KN=1,NEL
c$$$               CBEM(IN,JN)=CBEM(IN,JN)+BBEM(IN,KN)*EBEM(KN,JN)
c$$$            ENDDO
c$$$         ENDDO
c$$$      ENDDO
 
      IF(ILVOUT>=3) THEN
        WRITE (*,'(//,a)') ' ebem matrix = bbem*cbem-1'
        IF(NEL < 11) THEN
          DO IN = 1,NEL
             WRITE (*,'(10e13.5)')  (CBEM(IN,JN),JN=1,NEL)
          ENDDO
        ELSE
          DO IN = 1,10
             WRITE (*,'(10e13.5)')  (CBEM(IN,JN),JN=1,10)
          ENDDO
          WRITE (*,'(//,a)') ' ebem matrix e1i'
          WRITE (*,'(10e13.5)')  (CBEM(1,JN),JN=1,NEL)
        ENDIF
      ENDIF
 
      DO IN=1,NEL
         DO JN=1,NEL
            MFLE(IN,JN)=HALF*RHO*AF(IN)*(CBEM(IN,JN)+CBEM(JN,IN))
            BBEM(IN,JN)=MFLE(IN,JN)
         ENDDO
      ENDDO    
C
      IF(ILVOUT>=3) THEN
        WRITE (*,'(//,a)') ' fluid mass matrix'
        IF(NEL < 11) THEN
          DO IN = 1,NEL
             WRITE (*,'(10e13.5)')  (MFLE(IN,JN),JN=1,NEL)
          END DO
        ELSE
          DO IN = 1,10
             WRITE (*,'(10e13.5)')  (MFLE(IN,JN),JN=1,10)
          END DO
          WRITE (*,'(//,a)') ' fluid mass matrix m1i'
          WRITE (*,'(10e13.5)')  (MFLE(1,JN),JN=1,NEL)
        ENDIF
      ENDIF
C---------------------------------
C     COMPUTE RIGD BODY FLUID MASS
C---------------------------------
      MASSX=ZERO
      MASSY=ZERO
      MASSZ=ZERO
      DO IN=1,NEL
         DO JN=1,NEL
            MASSX=MASSX+NORMAL(1,IN)*MFLE(IN,JN)*NORMAL(1,JN)
            MASSY=MASSY+NORMAL(2,IN)*MFLE(IN,JN)*NORMAL(2,JN)
            MASSZ=MASSZ+NORMAL(3,IN)*MFLE(IN,JN)*NORMAL(3,JN)
         ENDDO    
      ENDDO    
      WRITE (IOUT,'(/7x,a,e13.5)')  'daa : rigid body fluid mass xx', MASSX
      WRITE (IOUT,'( 7x,a,e13.5)')  'daa : rigid body fluid mass yy', MASSY
      WRITE (IOUT,'( 7x,a,e13.5)')  'daa : rigid body fluid mass zz', MASSZ
 
C---------------------------------------------------
C 3- Compute inverse of fluid mass matrix
C---------------------------------------------------
      IF(KFORM==1) THEN
        DO IN=1,NEL
           DO JN=1,NEL
              EBEM(IN,JN)=ZERO
           ENDDO
        ENDDO
        CALL INVERT(BBEM,EBEM,NEL,ERR)
 
C
        IF(ILVOUT>=3) THEN
          WRITE (*,'(//,a)') ' inverse fluid mass matrix'
          IF(NEL < 11) THEN
            DO IN = 1,NEL
               WRITE (*,'(10e13.5)')  (EBEM(IN,JN),JN=1,NEL)
            END DO
          ELSE
            DO IN = 1,10
               WRITE (*,'(10e13.5)')  (EBEM(IN,JN),JN=1,10)
            END DO
            WRITE (*,'(//,a)') ' inverse fluid mass matrix mii'
            WRITE (*,'(10e13.5)')  (EBEM(JN,JN),JN=1,NEL)
          ENDIF
        ENDIF
 
        DO IN=1,NEL
           DO JN=1,NEL
              MFLE(IN,JN)=EBEM(IN,JN)
           ENDDO
        ENDDO
      ENDIF
      IF (ALLOCATED(BBEM)) DEALLOCATE(BBEM)
      IF (ALLOCATED(EBEM))DEALLOCATE(EBEM)    
            
      RETURN
      END