diri04.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE diri04
!                   *****************
!
     &(x1,x2,a11,a12,a21,a22,sm1,sm2,t1,t2,t3,t4,
     & xbor1,xbor2,lidir1,lidir2,mesh,kdir,msk,maskpt)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    TREATS THE DIRICHLET POINTS FOR THE FOLLOWING
!+                SYSTEM (BLOCK OF 4 MATRICES):
!code
!+         (     A11          A12              )  ( X1 )   ( SM1 )
!+         (                                   )  (    ) = (     )
!+         (     A21          A22              )  ( X2 )   ( SM2 )
!
!history  J-M HERVOUET (LNH)
!+        30/01/95
!+        V5P1
!+
!
!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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A12            |<->| MATRIX IN THE 2x2 LINEAR SYSTEM
!| A12            |<->| MATRIX IN THE 2x2 LINEAR SYSTEM
!| A21            |<->| MATRIX IN THE 2x2 LINEAR SYSTEM
!| A22            |<->| MATRIX IN THE 2x2 LINEAR SYSTEM
!| KDIR           |-->| CONVENTION FOR DIRICHLET BOUNDARY CONDITIONS
!| LIDIR1         |-->| TYPES OF BOUNDARY CONDITIONS FOR VARIABLE 1
!|                |   | IF LIMDIR(K) = KDIR THE KTH BOUNDARY POINT
!|                |   | IS OF DIRICHLET TYPE.
!| LIDIR2         |-->| TYPES OF BOUNDARY CONDITIONS FOR VARIABLE 2
!|                |   | IF LIMDIR(K) = KDIR THE 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.
!| SM1            |-->| FIRST RIGHT-HAND SIDE OF THE SYSTEM.
!| SM2            |-->| SECOND RIGHT-HAND SIDE OF THE SYSTEM.
!| T1             |<->| WORK DOUBLE PRECISION ARRAY IN A BIEF_OBJ
!| T2             |<->| WORK DOUBLE PRECISION ARRAY IN A BIEF_OBJ
!| T3             |<->| WORK DOUBLE PRECISION ARRAY IN A BIEF_OBJ
!| T4             |<->| WORK DOUBLE PRECISION ARRAY IN A BIEF_OBJ
!| XBOR1          |-->| DIRICHLET BOUNDARY CONDITIONS OF VARIABLE 1
!| XBOR2          |-->| DIRICHLET BOUNDARY CONDITIONS OF VARIABLE 2
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_diri04 => diri04
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj), INTENT(INOUT) :: X1,X2,SM1,SM2,T1,T2,T3,T4
      TYPE(bief_obj), INTENT(INOUT) :: A11,A12,A21,A22
      TYPE(bief_obj), INTENT(IN)    :: XBOR1,XBOR2,MASKPT
      INTEGER, INTENT(IN)           :: KDIR,LIDIR1(*),LIDIR2(*)
      TYPE(bief_mesh), INTENT(INOUT):: MESH
      LOGICAL, INTENT(IN)           :: MSK
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION C
!
      CHARACTER(LEN=1) STODIA
!
!-----------------------------------------------------------------------
!
! 1) BUILDS ARRAYS T1,T2 CONTAINING:
!    THE X1 AND X2 IMPOSED VALUES IF THE POINT IS OF TYPE DIRICHLET
!    0 OTHERWISE
!
!    X1,X2  3 TAKE THEIR DIRICHLET VALUE
!
!=======================================================================
!
!   BOUNDARY CONDITION FOR X1 : "XBOR1" IMPOSED
!
      CALL cpstvc(x1,t1)
      CALL os('X=C     ', x=t1, c=0.d0)
      CALL osdbif ( 'X=Y     ',t1,xbor1,lidir1,kdir,mesh)
!
!-----------------------------------------------------------------------
!
!   BOUNDARY CONDITIONS FOR X2 : "XBOR2" IMPOSED
!
      CALL cpstvc(x2,t2)
      CALL os  ('X=C     ', x=t2, c=0.d0)
      CALL osdbif ( 'X=Y     ',t2,xbor2,lidir2,kdir,mesh)
!
!=======================================================================
!
!   2) COMPUTES THE PRODUCT OF THE MATRIX FOR THE SYSTEM TO SOLVE
!      AND T1,T2
!      THE RESULT IS DEDUCTED FROM THE SECOND MEMBERS
!
      CALL matvec('X=AY    ',t3,a11,t1,c,mesh,lego=.false.)
      CALL matvec('X=X+AY  ',t3,a12,t2,c,mesh,lego=.true. )
      CALL matvec('X=AY    ',t4,a21,t1,c,mesh,lego=.false.)
      CALL matvec('X=X+AY  ',t4,a22,t2,c,mesh,lego=.true. )
!
      CALL cpstvc(x1,sm1)
      CALL cpstvc(x2,sm2)
      CALL os('X=X-Y   ', x=sm1, y=t3)
      CALL os('X=X-Y   ', x=sm2, y=t4)
!
!=======================================================================
!
!  SECOND MEMBERS OF THE EQUATIONS FOR DIRICHLET POINTS
!  PREPARES THE LINEAR SYSTEM
!
      CALL diraux(sm1,a11%D,xbor1,t1,x1,lidir1,kdir,mesh )
      CALL diraux(sm2,a22%D,xbor2,t2,x2,lidir2,kdir,mesh )
!
      IF(msk) THEN
        CALL ov('X=XY    ', x=sm1%R, y=maskpt%R, dim1=sm1%DIM1)
        CALL ov('X=XY    ', x=x1%R, y=maskpt%R, dim1=x1%DIM1)
        CALL ov('X=XY    ', x=t1%R, y=maskpt%R, dim1=t1%DIM1)
        CALL ov('X=XY    ', x=sm2%R, y=maskpt%R, dim1=sm2%DIM1)
        CALL ov('X=XY    ', x=x2%R, y=maskpt%R, dim1=x2%DIM1)
        CALL ov('X=XY    ', x=t2%R, y=maskpt%R, dim1=t2%DIM1)
      ENDIF
!
!=======================================================================
!
!   ERASES THE LINES AND COLUMNS FOR DIRICHLET POINTS
!
!   IT'S EQUIVALENT TO A DIAGONAL PRECONDITIONING WITH ARRAYS
!   T1,T2,T3
!
!   DOES NOT ALTER A11,A22,A33 DIAGONALS
!   BY GIVING THEM A DUMMY TYPE : '0'
!
!
!=======================================================================
! A11 PRECONDITIONING :
!=======================================================================
!
      stodia = a11%TYPDIA
      a11%TYPDIA='0'
      CALL om('M=DMD   ', m=a11, d=t1, mesh=mesh)
      a11%TYPDIA=stodia
!
!=======================================================================
! A12 PRECONDITIONING :
!=======================================================================
!
      CALL om('M=DM    ', m=a12, d=t1, mesh=mesh)
      CALL om('M=MD    ', m=a12, d=t2, mesh=mesh)
!
!=======================================================================
! A21 PRECONDITIONING :
!=======================================================================
!
      CALL om('M=DM    ', m=a21, d=t2, mesh=mesh)
      CALL om('M=MD    ', m=a21, d=t1, mesh=mesh)
!
!=======================================================================
! A22 PRECONDITIONING :
!=======================================================================
!
      stodia = a22%TYPDIA
      a22%TYPDIA='0'
      CALL om('M=DMD   ', m=a22, d=t2, mesh=mesh)
      a22%TYPDIA=stodia
!
!-----------------------------------------------------------------------
!
      RETURN
      END