preverebe.f

Go to the documentation of this file.

!                   ********************
                    SUBROUTINE preverebe
!                   ********************
!
     &(xaux,ax,typdia,typext,ikle,npoin,nelem,nelmax,mesh,typemesh)
!
!***********************************************************************
! BIEF   V6P2                                   21/08/2010
!***********************************************************************
!
!brief    BUILDS TRIDIAGONAL SYSTEMS FOR EVERY VERTICAL,
!+                BY LUMPING A MATRIX DEFINED ON PRISMS.
!
!warning  THE JACOBIAN MUST BE POSITIVE
!
!history  J-M HERVOUET (LNHE)
!+        02/06/08
!+        V5P9
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AX             |-->| MATRIX OFF-DIAGONAL TERMS
!| IKLE           |-->| CONNECTIVITY TABLE.
!| MESH           |-->| MESH STRUCTURE
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| TYPDIA         |-->| TYPE OF DIAGONAL:
!|                |   | TYPDIA = 'Q' : ANY VALUE
!|                |   | TYPDIA = 'I' : IDENTITY
!|                |   | TYPDIA = '0' : ZERO
!| TYPEXT         |-->| TYPE OF OFF-DIAGONAL TERMS
!|                |   | TYPEXT = 'Q' : ANY VALUE
!|                |   | TYPEXT = 'S' : SYMMETRIC
!|                |   | TYPEXT = '0' : ZERO
!| TYPEMESH       |-->| TYPE OF MESH (40: PRISMS, 50: PRISMS CUT INTO
!|                |   | TETRAHEDRONS)
!| XAUX           |<--| TRIDIAGONAL MATRIX
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_preverebe => preverebe
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX,NPOIN,TYPEMESH
      INTEGER, INTENT(IN) :: IKLE(nelmax,6)
!
      DOUBLE PRECISION, INTENT(IN) :: AX(nelmax,*)
      DOUBLE PRECISION, INTENT(INOUT) :: XAUX(npoin,*)
!
      CHARACTER(LEN=1), INTENT(IN) :: TYPDIA,TYPEXT
!
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I1,I2,I3,I4,I5,I6,IELEM,NPLAN,IAN,ICOM,NPOIN2
!
!-----------------------------------------------------------------------
!
!     HERE WE CONSIDER THAT NPOIN < NELMAX TO USE XAUX AS XAUX(NPOIN,3)
!
!     XAUX(I,1) IS COEFFICIENT OF POINT BELOW I IN EQUATION OF POINT I
!     XAUX(I,2) IS THE DIAGONAL
!     XAUX(I,3) IS COEFFICIENT OF POINT ABOVE I IN EQUATION OF POINT I
!
!-----------------------------------------------------------------------
!     INITIALISES THE DIAGONAL AND OFF-DIAGONAL TERMS
!-----------------------------------------------------------------------
!
!     OFF-DIAGONAL
!
      CALL ov('X=C     ', x=xaux(1,1), c=0.d0, dim1=npoin)
      CALL ov('X=C     ', xaux(1,3), c=0.d0, dim1=npoin)
!
!     DIAGONAL
!
      IF(typdia(1:1).EQ.'0') THEN
        CALL ov('X=C     ', x=xaux(1,2), c=0.d0, dim1=npoin)
      ELSEIF(typdia(1:1).EQ.'I') THEN
        CALL ov('X=C     ', x=xaux(1,2), c=1.d0, dim1=npoin)
      ELSEIF(typdia(1:1).EQ.'Q') THEN
        CALL ov('X=Y     ', x=xaux(1,2), c=0.d0, dim1=npoin)
      ELSE
        WRITE(lu,*) typdia
        WRITE(lu,*) 'UNKNOWN TYPE OF DIAGONAL IN PREVEREBE'
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!     LUMPS THE OFF-DIAGONAL TERMS
!
!     ONLY VERTICAL OF POINT, LUMPING ALL TERMS HAS BEEN TESTED, SEE
!     LINE COMMENTED, BUT DOES NOT WORK (E.G. TRY SOLITARY WAVE)
!
!
!-----------------------------------------------------------------------
!
      IF(typemesh.EQ.40) THEN
!
      IF(typext.EQ.'Q') THEN
        DO ielem=1,nelem
          i1=ikle(ielem,1)
          i2=ikle(ielem,2)
          i3=ikle(ielem,3)
          i4=ikle(ielem,4)
          i5=ikle(ielem,5)
          i6=ikle(ielem,6)
          xaux(i1,3)=xaux(i1,3)+ax(ielem,03) ! TERM 1-4
          xaux(i2,3)=xaux(i2,3)+ax(ielem,08) ! TERM 2-5
          xaux(i3,3)=xaux(i3,3)+ax(ielem,12) ! TERM 3-6
          xaux(i4,1)=xaux(i4,1)+ax(ielem,18) ! TERM 4-1
          xaux(i5,1)=xaux(i5,1)+ax(ielem,23) ! TERM 5-2
          xaux(i6,1)=xaux(i6,1)+ax(ielem,27) ! TERM 6-3
        ENDDO
      ELSEIF(typext.EQ.'S') THEN
        DO ielem=1,nelem
          i1=ikle(ielem,1)
          i2=ikle(ielem,2)
          i3=ikle(ielem,3)
          i4=ikle(ielem,4)
          i5=ikle(ielem,5)
          i6=ikle(ielem,6)
!         XAUX(I1,2)=XAUX(I1,2)+AX(IELEM,01) ! TERM 1-2
!         XAUX(I2,2)=XAUX(I2,2)+AX(IELEM,01) ! TERM 2-1
!         XAUX(I1,2)=XAUX(I1,2)+AX(IELEM,02) ! TERM 1-3
!         XAUX(I2,2)=XAUX(I2,2)+AX(IELEM,02) ! TERM 3-1
          xaux(i1,3)=xaux(i1,3)+ax(ielem,03) ! TERM 1-4
          xaux(i4,1)=xaux(i4,1)+ax(ielem,03) ! TERM 4-1
!         XAUX(I1,3)=XAUX(I1,3)+AX(IELEM,04) ! TERM 1-5
!         XAUX(I5,1)=XAUX(I5,1)+AX(IELEM,04) ! TERM 5-1
!         XAUX(I1,3)=XAUX(I1,3)+AX(IELEM,05) ! TERM 1-6
!         XAUX(I6,1)=XAUX(I6,1)+AX(IELEM,05) ! TERM 6-1
!         XAUX(I2,2)=XAUX(I2,2)+AX(IELEM,06) ! TERM 2-3
!         XAUX(I3,2)=XAUX(I3,2)+AX(IELEM,06) ! TERM 3-2
!         XAUX(I2,3)=XAUX(I2,3)+AX(IELEM,07) ! TERM 2-4
!         XAUX(I4,1)=XAUX(I4,1)+AX(IELEM,07) ! TERM 4-2
          xaux(i2,3)=xaux(i2,3)+ax(ielem,08) ! TERM 2-5
          xaux(i5,1)=xaux(i5,1)+ax(ielem,08) ! TERM 5-2
!         XAUX(I2,3)=XAUX(I2,3)+AX(IELEM,09) ! TERM 2-6
!         XAUX(I6,1)=XAUX(I6,1)+AX(IELEM,09) ! TERM 6-2
!         XAUX(I3,3)=XAUX(I3,3)+AX(IELEM,10) ! TERM 3-4
!         XAUX(I4,1)=XAUX(I4,1)+AX(IELEM,10) ! TERM 4-3
!         XAUX(I3,3)=XAUX(I3,3)+AX(IELEM,11) ! TERM 3-5
!         XAUX(I5,1)=XAUX(I5,1)+AX(IELEM,11) ! TERM 5-3
          xaux(i3,3)=xaux(i3,3)+ax(ielem,12) ! TERM 3-6
          xaux(i6,1)=xaux(i6,1)+ax(ielem,12) ! TERM 6-3
!         XAUX(I4,2)=XAUX(I4,2)+AX(IELEM,13) ! TERM 4-5
!         XAUX(I5,2)=XAUX(I5,2)+AX(IELEM,13) ! TERM 5-4
!         XAUX(I4,2)=XAUX(I4,2)+AX(IELEM,14) ! TERM 4-6
!         XAUX(I6,2)=XAUX(I6,2)+AX(IELEM,14) ! TERM 6-4
!         XAUX(I5,2)=XAUX(I5,2)+AX(IELEM,15) ! TERM 5-6
!         XAUX(I6,2)=XAUX(I6,2)+AX(IELEM,15) ! TERM 6-5
        ENDDO
      ELSEIF(typext.EQ.'0') THEN
!       NOTHING TO DO (BUT WHAT'S THE USE OF AN ITERATIVE SOLVER ?)
      ELSE
        WRITE(lu,*) typext
        WRITE(lu,*) 'UNKNOWN TYPE OF OFF-DIAGONAL TERMS'
        WRITE(lu,*) 'IN PREVEREBE'
        CALL plante(1)
        stop
      ENDIF
!
      ELSE
        WRITE(lu,*) typemesh
        WRITE(lu,*) 'UNKNOWN TYPE OF MESH'
        WRITE(lu,*) 'IN PREVEREBE: ',typemesh
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!     PARALLEL MODE
!
      IF(ncsize.GT.1) THEN
        ian    = 3
        icom   = 2
        npoin2 = bief_nbpts(11,mesh)
        nplan=npoin/npoin2
        CALL parcom2(xaux(1,1),xaux(1,2),xaux(1,3),
     &               npoin2,nplan,icom,ian,mesh)
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END