descen.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE descen
!                   *****************
!
     &(x, xa,typexa,b,ikle,nelem,nelmax,npoin,ielm,ditr,copy,lv)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    SOLVES THE SYSTEM L X = B (ELEMENT BY ELEMENT).
!code
!+            MATRIX L IS HERE THE RESULT OF THE FACTORISATION
!+            PERFORMED IN SUBROUTINE DECLDU.
!+
!+            EACH ELEMENTARY MATRIX IS DECOMPOSED IN THE FORM:
!+
!+            LE X DE X UE
!+
!+            LE : LOWER TRIANGULAR WITH 1S ON THE DIAGONAL
!+            DE : DIAGONAL
!+            UE : UPPER TRIANGULAR WITH 1S ON THE DIAGONAL
!+
!+                                                T
!+            IF THE MATRIX IS SYMMETRICAL : LE =  UE
!
!code
!+-----------------------------------------------------------------------
!+  MEANING OF IELM :
!+
!+  TYPE OF ELEMENT      NUMBER OF POINTS      CODED IN THIS SUBROUTINE
!+
!+  11 : P1 TRIANGLE            3                       YES
!+  12 : P2 TRIANGLE            6
!+  13 : P1-ISO P1 TRIANGLE     6
!+  14 : P2 TRIANGLE            7
!+  21 : Q1 QUADRILATERAL       4                       YES
!+  22 : Q2 QUADRILATERAL       8
!+  24 : Q2 QUADRILATERAL       9
!+  31 : P1 TETRAHEDRON         4
!+  32 : P2 TETRAHEDRON        10
!+  41 : TELEMAC-3D PRISMS      6                       YES
!+  41 : PRISMS SPLIT IN TETRAHEDRONS     4             YES
!+
!+-----------------------------------------------------------------------
!
!history  J-M HERVOUET (LNH)    ; F  LEPEINTRE (LNH)
!+        18/08/94
!+        V5P3
!+
!
!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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| B              |<--| RIGHT-HAND SIDE OF THE LINEAR SYSTEM TO BE SOLVED
!| COPY           |-->| IF .TRUE. B IS COPIED INTO X TO START WITH
!| DITR           |-->| CHARACTER, IF  'D' : DIRECT MATRIX A CONSIDERED
!|                |   |                'T' : TRANSPOSED MATRIX A CONSIDERED
!| IELM           |-->| TYPE OF ELEMENT (SEE ABOVE)
!| IKLE           |-->| CONNECTIVITY TABLE.
!| LV             |-->| VECTOR LENGTH OF THE MACHINE
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| TYPEXA         |-->| TYPE OF OFF-DIAGONAL TERMS IN THE MATRIX
!| X              |<--| SOLUTION OF THE SYSTEM AX = B
!| XA             |<--| OFF-DIAGONAL TERMS OF THE MATRIX
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_descen => descen
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: IELM,NPOIN,NELEM,NELMAX,LV
      INTEGER         , INTENT(IN)    :: IKLE(nelmax,*)
      DOUBLE PRECISION, INTENT(INOUT) :: X(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: XA(nelmax,*),B(npoin)
      CHARACTER(LEN=1), INTENT(IN)    :: TYPEXA,DITR
      LOGICAL         , INTENT(IN)    :: COPY
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
! 1) INITIALISES : X = SECOND MEMBER
!
      IF(copy) CALL ov('X=Y     ', x=x, y=b, dim1=npoin)
!
!-----------------------------------------------------------------------
!
! 2) INVERTS THE LOWER TRIANGULAR MATRICES (DESCENT)
!
!     2.1) TRANSPOSE CASE
!
      IF(typexa(1:1).EQ.'S' .OR.
     &  (typexa(1:1).EQ.'Q'.AND.ditr(1:1).EQ.'T')) THEN
!
      IF(ielm.EQ.11) THEN
!
        CALL des11(x,xa(1,1),xa(1,2),xa(1,3),
     &             ikle(1,1),ikle(1,2),ikle(1,3),nelem,nelmax,npoin,lv)
!
      ELSEIF(ielm.EQ.21.OR.ielm.EQ.12.OR.ielm.EQ.31.OR.ielm.EQ.51) THEN
!
        CALL des21(x,xa(1,1),xa(1,2),xa(1,3),xa(1,4),xa(1,5),xa(1,6),
     &             ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &             nelem,nelmax,npoin,lv)
!
      ELSEIF(ielm.EQ.41) THEN
!
        CALL des41(x,xa(1,1),xa(1,2),xa(1,3),xa(1,4),xa(1,5),xa(1,6),
     &             xa(1,7),xa(1,8),xa(1,9),xa(1,10),xa(1,11),xa(1,12),
     &             xa(1,13),xa(1,14),xa(1,15),
     &             ikle(1,1),ikle(1,2),ikle(1,3),
     &             ikle(1,4),ikle(1,5),ikle(1,6),nelem,nelmax,npoin,lv)
!
!  IELM NOT IMPLEMENTED: ERROR
!
      ELSE
!
        WRITE(lu,101) ielm
101   FORMAT(1x,'DESCEN (BIEF) : IELM = ',1i6,' ELEMENT NOT AVAILABLE')
        CALL plante(1)
        stop
!
      ENDIF
!
!     2.2) DIRECT CASE
!
      ELSEIF(typexa(1:1).EQ.'Q'.AND.ditr(1:1).EQ.'D') THEN
!
      IF(ielm.EQ.11) THEN
!
        CALL des11(x,xa(1,4),xa(1,5),xa(1,6),
     &             ikle(1,1),ikle(1,2),ikle(1,3),nelem,nelmax,npoin,lv)
!
      ELSEIF(ielm.EQ.21.OR.ielm.EQ.12.OR.ielm.EQ.31.OR.ielm.EQ.51) THEN
!
        CALL des21(x,xa(1,7),xa(1,8),xa(1,9),xa(1,10),xa(1,11),xa(1,12),
     &             ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &             nelem,nelmax,npoin,lv)
!
      ELSEIF(ielm.EQ.41) THEN
!
        CALL des41(x,
     &            xa(1,16),xa(1,17),xa(1,18),xa(1,19),xa(1,20),xa(1,21),
     &            xa(1,22),xa(1,23),xa(1,24),xa(1,25),xa(1,26),xa(1,27),
     &            xa(1,28),xa(1,29),xa(1,30),
     &            ikle(1,1),ikle(1,2),ikle(1,3),
     &            ikle(1,4),ikle(1,5),ikle(1,6),nelem,nelmax,npoin,lv)
!
!  IELM NOT IMPLEMENTED: ERROR
!
      ELSE
!
        WRITE(lu,101) ielm
        CALL plante(1)
        stop
!
      ENDIF
!
!     2.3) CASE NOT IMPLEMENTED
!
      ELSE
        WRITE(lu,201) typexa(1:1),ditr(1:1)
201     FORMAT(1x,'DESCEN (BIEF) : UNEXPECTED TYPE OF MATRIX :',a1,/,
     &         1x,'WITH DITR=',a1)
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END