diri01.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE diri01
!                   *****************
!
     &(f, s, sm ,fbor,limdir,work1,work2,mesh,kdir,msk,maskpt)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    TAKES INTO ACCOUNT POINTS OF TYPE DIRICHLET IN A SYSTEM
!+                OF LINEAR EQUATIONS WITH SYMMETRICAL MATRIX.
!+
!+            IN THE EQUATIONS FOR POINTS NOT OF TYPE DIRICHLET :
!+                DIRICHLET VALUES ARE REMOVED.
!+
!+            IN THE EQUATIONS FOR POINTS OF TYPE DIRICHLET :
!+                DEFINES AN EQUATION FIXING THE IMPOSED VALUE.
!
!warning  THIS SUBROUTINE IS NOT PROTECTED AGAINST DIAGONAL EQUAL
!+            TO 0 ON DIRICHLET POINTS; IT WILL THEN SET AN EQUATION
!+            0 X = 0 ON SUCH POINTS
!
!history  J-M HERVOUET (LNHE)
!+        07/08/2009
!+        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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| F              |<->| VALUES AT TIME N+1 AND INITIALISATION
!| FBOR           |-->| DIRICHLET BOUNDARY CONDITIONS
!| KDIR           |-->| CONVENTION FOR DIRICHLET BOUNDARY CONDITIONS
!| LIMDIR         |-->| TYPES OF BOUNDARY CONDITIONS
!|                |   | IF LIMDIR(K) = KDIR LE KTH BOUNDARY POINT
!|                |   | IS OF DIRICHLET TYPE.
!| MASKPT         |-->| MASKING PER POINT.
!|                |   | =1. : NORMAL   =0. : MASKED
!| MESH           |-->| MESH STRUCTURE
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
!| S              |<->| MATRIX OF THE LINEAR SYSTEM
!| SM             |-->| RIGHT HAND SIDE
!| WORK1          |<->| WORK ARRAY
!| WORK2          |<->| WORK ARRAY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_diri01 => diri01
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj), INTENT(INOUT) :: F,S,SM,WORK1,WORK2
      TYPE(bief_obj), INTENT(IN)    :: FBOR,MASKPT
      INTEGER, INTENT(IN) :: LIMDIR(*), KDIR
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
      LOGICAL, INTENT(IN) :: MSK
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION C
!
      INTEGER IELMSM,IELMFB
!
      CHARACTER(LEN=1) OLDDIA
!
!----------------------------------------------------------------------
!
!     DEPLOYS THE MESH STRUCTURE
!
!----------------------------------------------------------------------
!
!     BUILDS AN ARRAY WITH 0S EVERYWHERE EXCEPT AT DIRICHLET POINTS
!     FOR WHICH THE VALUE IS TAKEN FROM FBOR
!     FBOR MUST BE 0 WHEN THE POINT IS NOT OF TYPE DIRICHLET
!
      CALL cpstvc(sm,work1)
!
      ielmsm=sm%ELM
      ielmfb=fbor%ELM
      IF(ielmsm.EQ.ielmfb) THEN
        CALL matvec( 'X=AY    ' ,work2,s,fbor,c, mesh )
      ELSE
        CALL os( 'X=0     ' , x=work1 )
        CALL osdbif( 'X=Y     ' ,work1,fbor,limdir,kdir,mesh)
        CALL matvec( 'X=AY    ' ,work2,s,work1,c, mesh )
      ENDIF
!
!----------------------------------------------------------------------
!
!     THE PRODUCT S WORK1 IS DEDUCTED FROM THE SECOND MEMBER.
!     IT MEANS THAT THE VALUES AT DIRICHLET POINTS ARE NO LONGER
!     UNKNOWNS IN THE EQUATIONS FOR THE OTHER POINTS.
!
      CALL os( 'X=X-Y   ' , x=sm , y=work2 )
!
!----------------------------------------------------------------------
!
!     BUILDS AN ARRAY WITH 1S EVERYWHERE EXCEPT AT DIRICHLET POINTS
!     FOR WHICH IT'S 0
!
!     WHAT'S MORE, AN EQUATION OF THE FORM DS(N) * X = DS(N) * FBOR
!     (WILL GIVE X=FBOR) IS SET IN THE MATRIX FOR DIRICHLET POINTS;
!     AND F IS INITIALISED TO ITS KNOWN VALUE.
!     THIS ASSUMES THAT DS(N) IS NOT 0
!
      CALL diraux(sm,s%D,fbor,work2,f,limdir,kdir,mesh )
!
!     MASKING : FOR THE POINTS OF MASKED ELEMENTS THE EQUATION X=0
!               IS SET FOR THE DIAGONAL COEFFICIENT PRES
!
      IF(msk) THEN
        CALL ov('X=XY    ', x=sm%R   , y=maskpt%R, dim1=sm%DIM1)
        CALL ov('X=XY    ', x=f%R    , y=maskpt%R, dim1=f%DIM1)
        CALL ov('X=XY    ', x=work2%R, y=maskpt%R, dim1=work2%DIM1)
      ENDIF
!
!----------------------------------------------------------------------
!
!     WORK2 * S * WORK2 :
!     ERASES THE LINES AND COLUMNS IN S WHICH CORRESPOND TO DIRICHLET
!     POINTS
!     DOES NOT ALTER THE DIAGONAL BY DECLARING IT 0 HERE
!
      olddia=s%TYPDIA
      s%TYPDIA='0'
      CALL om('M=DMD   ', m=s, d=work2, mesh=mesh)
      s%TYPDIA=olddia
!
!----------------------------------------------------------------------
!
      RETURN
      END