preverseg.f

Go to the documentation of this file.

!                   ********************
                    SUBROUTINE preverseg
!                   ********************
!
     &(xaux,ad,ax,typdia,typext,npoin,mesh,nseg3d,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  JMH
!+        11/08/09
!+
!+   CROSSED AND VERTICAL SEGMENTS SWAPPED (SEE STOSEG41)
!
!history  J-M HERVOUET (LNHE)
!+        11/08/09
!+        V6P0
!+
!
!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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AD             |-->| MATRIX DIAGONAL
!| AX             |-->| MATRIX OFF-DIAGONAL TERMS
!| MESH           |-->| MESH STRUCTURE
!| NPOIN          |-->| NUMBER OF POINTS
!| NSEG3D         |-->| NUMBER OF SEGMENTS IN 3D MESH
!| 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_preverseg => preverseg
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NPOIN,NSEG3D,TYPEMESH
!
      DOUBLE PRECISION, INTENT(IN)    :: AD(npoin),AX(nseg3d,*)
      DOUBLE PRECISION, INTENT(INOUT) :: XAUX(npoin,*)
!
      CHARACTER(LEN=1), INTENT(IN) :: TYPDIA,TYPEXT
!
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I2,I3,NPLAN,IAN,ICOM,NPOIN2,SEGUP,SEGDOWN,NSEG2D
      INTEGER IPLAN,NSEGH
!
!-----------------------------------------------------------------------
!
!     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
!-----------------------------------------------------------------------
!
      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), y=ad, dim1=npoin)
      ELSE
        WRITE(lu,*) typdia
        WRITE(lu,*) 'UNKNOWN TYPE OF DIAGONAL IN PREVERSEG'
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!     LUMPS THE OFF-DIAGONAL TERMS CORRESPONDING TO VERTICAL SEGMENTS
!-----------------------------------------------------------------------
!
      npoin2 = bief_nbpts(11,mesh)
      nplan  = npoin/npoin2
      nseg2d = bief_nbseg(11,mesh)
      nsegh  = nseg2d*nplan
!
      IF(typemesh.EQ.40.OR.typemesh.EQ.50) THEN
!
      IF(typext.EQ.'Q') THEN
!       PLANE ON THE BOTTOM
        DO i2=1,npoin2
          segup=nsegh+i2
          xaux(i2,1)=0.d0
          xaux(i2,3)=ax(segup,1)
        ENDDO
!       PLANE AT THE FREE SURFACE
        DO i2=1,npoin2
          i3=i2+(nplan-1)*npoin2
          segdown=nsegh+npoin2*(nplan-2)+i2
          xaux(i3,1)=ax(segdown,2)
          xaux(i3,3)=0.d0
        ENDDO
!       OTHER PLANES
        IF(nplan.GT.2) THEN
        DO iplan=2,nplan-1
          DO i2=1,npoin2
            i3=i2+(iplan-1)*npoin2
            segdown=nsegh+npoin2*(iplan-2)+i2
            segup  =segdown+npoin2
            xaux(i3,1)=ax(segdown,2)
            xaux(i3,3)=ax(segup,1)
          ENDDO
        ENDDO
        ENDIF
      ELSEIF(typext.EQ.'S') THEN
!       PLANE ON THE BOTTOM
        DO i2=1,npoin2
          segup=nsegh+i2
          xaux(i2,1)=0.d0
          xaux(i2,3)=ax(segup,1)
        ENDDO
!       PLANE AT THE FREE SURFACE
        DO i2=1,npoin2
          i3=i2+(nplan-1)*npoin2
          segdown=nsegh+npoin2*(nplan-2)+i2
          xaux(i3,1)=ax(segdown,1)
          xaux(i3,3)=0.d0
        ENDDO
!       OTHER PLANES
        IF(nplan.GT.2) THEN
        DO iplan=2,nplan-1
          DO i2=1,npoin2
            i3=i2+(iplan-1)*npoin2
            segdown=nsegh+npoin2*(iplan-2)+i2
            segup  =segdown+npoin2
            xaux(i3,1)=ax(segdown,1)
            xaux(i3,3)=ax(segup,1)
          ENDDO
        ENDDO
        ENDIF
      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 PREVERSEG'
        CALL plante(1)
        stop
      ENDIF
!
      ELSE
        WRITE(lu,*) typemesh
        WRITE(lu,*) 'UNKNOWN TYPE OF MESH'
        WRITE(lu,*) 'IN PREVERSEG: ',typemesh
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!     PARALLEL MODE
!
      IF(ncsize.GT.1) THEN
        ian    = 3
        icom   = 2
        CALL parcom2(xaux(1,1),xaux(1,2),xaux(1,3),
     &               npoin2,nplan,icom,ian,mesh)
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END