gettriseg.f

Go to the documentation of this file.

!                   ********************
                    SUBROUTINE gettriseg
!                   ********************
!
     &(xaux,ad,ax,teta,npoin,mesh,nseg3d,nseg2d,nplan,npoin2,ielm3)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    GETS THE TRIDIAGONAL PART OF A DIFFUSION MATRIX ON
!+                PRISMS AND REMOVES IT FROM THE INITIAL MATRIX.
!code
!+           IF MTRI IS THIS TRIDIAGONAL PART, MAUX THE RESULT AND MDIF
!+           THE DIFFUSION MATRIX, THIS SUBROUTINE DOES:
!+
!+           MAUX = TETA * MTRI
!+           MDIF CHANGED INTO (1-TETA) * MDIF
!+
!+           SEGMENT STORAGE FOR MDIFF HERE !!!!!!!!!!!!!!!!!!!!!!!!
!
!warning  THE JACOBIAN MUST BE POSITIVE
!
!history  J-M HERVOUET (LNHE)
!+        11/08/09
!+        V6P0
!+   CROSSED AND VERTICAL SEGMENTS SWAPPED (SEE STOSEG41)
!
!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             |-->| DIAGONAL TERMS OF MATRIX
!| AX             |-->| OFF-DIAGONAL TERMS OF MATRIX
!|                |   | HERE DIMENSION 1 BECAUSE SYMMETRY
!| IELM3          |-->| TYPE OF ELEMENT
!| MESH           |-->| MESH STRUCTURE
!| NPLAN          |-->| NUMBER OF PLANES
!| NPOIN          |-->| NUMBER OF POINTS
!| NPOIN2         |-->| NUMBER OF POINTS OF 2D MESH
!| NSEG2D         |-->| NUMBER OF SEGMENTS IN 2D
!| NSEG3D         |-->| NUMBER OF SEGMENTS IN 3D
!| TETA           |-->| COEFFICIENT USED IN THE RESULT
!| XAUX           |<--| THE RESULTING MATRIX
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_gettriseg => gettriseg
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NPOIN,NSEG3D,NSEG2D,NPLAN,NPOIN2,IELM3
!
      DOUBLE PRECISION, INTENT(IN)    :: TETA
      DOUBLE PRECISION, INTENT(INOUT) :: XAUX(npoin,*),AX(nseg3d)
      DOUBLE PRECISION, INTENT(INOUT) :: AD(npoin)
!
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I2,I3,IPLAN,IAN,ICOM,SEGUP,SEGDOWN,NSEGH,NSEGV
!
!-----------------------------------------------------------------------
!
!     CONSIDERS HERE 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 TERMS
!-----------------------------------------------------------------------
!
      CALL ov('X=CY    ', x=xaux(1,2), y=ad, c=teta, dim1=npoin)
      CALL ov('X=CX    ', x=ad, c=1.d0-teta, dim1=npoin)
!
!-----------------------------------------------------------------------
!     TRIDIAGONAL TERMS
!-----------------------------------------------------------------------
!
!     VERTICAL SEGMENTS HAVE THE SAME POSITION, SEE STOSEG41,STOSEG51
!
      IF(ielm3.EQ.41.OR.ielm3.EQ.51) THEN
!
        nsegh=nseg2d*nplan
        nsegv=npoin2*(nplan-1)
!
!       PLANE ON THE BOTTOM
!
        DO i2=1,npoin2
          segup=nsegh+i2
          xaux(i2,1)=0.d0
          xaux(i2,3)=teta*ax(segup)
        ENDDO
!
!       PLANE AT THE FREE SURFACE
!
        DO i2=1,npoin2
          i3=i2+(nplan-1)*npoin2
          segdown=nsegh+npoin2*(nplan-2)+i2
          xaux(i3,1)=teta*ax(segdown)
          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)=teta*ax(segdown)
              xaux(i3,3)=teta*ax(segup)
            ENDDO
          ENDDO
        ENDIF
!
        CALL ov('X=CX    ', x=ax(nsegh+1:nsegh+nsegv),
     &          c=1.d0-teta, dim1=nsegv)
!
      ELSE
        WRITE(lu,*) 'GETTRISEG: UNKNOWN ELEMENT:',ielm3
        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