om1201.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE om1201
!                   *****************
!
     &(op ,  dm,typdim,xm,typexm,   dn,typdin,xn,typexn,   c,
     & nulone,nelbor,nbor,nelmax,ndiag,nptfr,nelebx,neleb)
!
!***********************************************************************
! BIEF   V7P0                                         21/08/2010
!***********************************************************************
!
!brief    OPERATIONS ON MATRICES.
!code
!+   M: QUASI-BUBBLE TRIANGLE
!+   N: BOUNDARY MATRIX
!+   D: DIAGONAL MATRIX
!+   C: CONSTANT
!+
!+   OP IS A STRING OF 8 CHARACTERS, WHICH INDICATES THE OPERATION TO BE
!+   PERFORMED ON MATRICES M AND N, D AND C.
!+
!+   THE RESULT IS MATRIX M.
!+
!+      OP = 'M=M+N   '  : ADDS N TO M
!
!code
!+  CONVENTION FOR THE STORAGE OF EXTRA-DIAGONAL TERMS:
!+
!+      XM(     , 1)  ---->  M(1,2)
!+      XM(     , 2)  ---->  M(1,3)
!+      XM(     , 3)  ---->  M(1,4)
!+      XM(     , 4)  ---->  M(2,3)
!+      XM(     , 5)  ---->  M(2,4)
!+      XM(     , 6)  ---->  M(3,4)
!+      XM(     , 7)  ---->  M(2,1)
!+      XM(     , 8)  ---->  M(3,1)
!+      XM(     , 9)  ---->  M(4,1)
!+      XM(     ,10)  ---->  M(3,2)
!+      XM(     ,11)  ---->  M(4,2)
!+      XM(     ,12)  ---->  M(4,3)
!
!history  J-M HERVOUET (LNHE)
!+        23/06/2008
!+        V5P9
!+
!
!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
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        13/03/2014
!+        V7P0
!+   Now written to enable different numbering of boundary points and
!+   boundary segments.
!
!history  S.E.BOURBAN (HRW)
!+        21/03/2017
!+        V7P3
!+   Replacement of the DATA declarations by the PARAMETER associates
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| C              |-->| A GIVEN CONSTANT USED IN OPERATION OP
!| DM             |<->| DIAGONAL OF M
!| DN             |-->| DIAGONAL OF N
!| NBOR           |-->| GLOBAL NUMBER OF BOUNDARY POINTS
!| NDIAG          |-->| NUMBER OF TERMS IN THE DIAGONAL
!| NELBOR         |-->| FOR THE KTH BOUNDARY EDGE, GIVES THE CORRESPONDING
!|                |   | ELEMENT.
!| NELEB          |-->| NUMBER OF BOUNDARY ELEMENTS
!| NELEBX         |-->| MAXIMUM NUMBER OF BOUNDARY ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| NULONE         |-->| GOES WITH ARRAY NELBOR. NELBOR GIVES THE
!|                |   | ADJACENT ELEMENT, NULONE GIVES THE LOCAL
!|                |   | NUMBER OF THE FIRST NODE OF THE BOUNDARY EDGE
!|                |   | I.E. 1, 2 OR 3 FOR TRIANGLES.
!| OP             |-->| OPERATION TO BE DONE (SEE ABOVE)
!| TYPDIM         |<->| TYPE OF DIAGONAL OF M:
!|                |   | TYPDIM = 'Q' : ANY VALUE
!|                |   | TYPDIM = 'I' : IDENTITY
!|                |   | TYPDIM = '0' : ZERO
!| TYPDIN         |<->| TYPE OF DIAGONAL OF N:
!|                |   | TYPDIN = 'Q' : ANY VALUE
!|                |   | TYPDIN = 'I' : IDENTITY
!|                |   | TYPDIN = '0' : ZERO
!| TYPEXM         |-->| TYPE OF OFF-DIAGONAL TERMS OF M:
!|                |   | TYPEXM = 'Q' : ANY VALUE
!|                |   | TYPEXM = 'S' : SYMMETRIC
!|                |   | TYPEXM = '0' : ZERO
!| TYPEXN         |-->| TYPE OF OFF-DIAGONAL TERMS OF N:
!|                |   | TYPEXN = 'Q' : ANY VALUE
!|                |   | TYPEXN = 'S' : SYMMETRIC
!|                |   | TYPEXN = '0' : ZERO
!| XM             |-->| OFF-DIAGONAL TERMS OF M
!| XN             |-->| OFF-DIAGONAL TERMS OF N
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_om1201 => om1201
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NELMAX,NDIAG,NPTFR,NELEBX,NELEB
      CHARACTER(LEN=8), INTENT(IN)    :: OP
      INTEGER, INTENT(IN)             :: NULONE(nelebx),NELBOR(nelebx)
      INTEGER, INTENT(IN)             :: NBOR(*)
      DOUBLE PRECISION, INTENT(IN)    :: DN(*),XN(*)
      DOUBLE PRECISION, INTENT(INOUT) :: DM(*),XM(nelmax,*)
      CHARACTER(LEN=1), INTENT(INOUT) :: TYPDIM,TYPEXM,TYPDIN,TYPEXN
      DOUBLE PRECISION, INTENT(IN)    :: C
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER K,IEL
!
      DOUBLE PRECISION Z(1)
!
!-----------------------------------------------------------------------
!     BEWARE: ONLY WORKS FOR QUASI-BUBBLE
      INTEGER :: CORNSY(4,2)
      parameter( cornsy = reshape( (/
     &                        1,4,8,0,  7,10,2,0/), shape=(/ 4,2 /) ) )
      INTEGER :: CORSYM(4)
      parameter( corsym = (/ 1,4,2,0      /) )
!
!-----------------------------------------------------------------------
!
      IF(op(1:8).EQ.'M=M+N   ') THEN
!
        IF(typdim.EQ.'Q'.AND.typdim.EQ.'Q'.AND.ndiag.GE.nptfr) THEN
          CALL ovdb( 'X=X+Y   ' , dm , dn , z , c , nbor , nptfr )
        ELSE
          WRITE(lu,199) typdim(1:1),op(1:8),typdin(1:1)
199       FORMAT(1x,'OM1201 (BIEF) : TYPDIM = ',a1,' NOT IMPLEMENTED',
     &      /,1x,'FOR THE OPERATION : ',a8,' WITH TYPDIN = ',a1)
          CALL plante(1)
          stop
        ENDIF
!
        IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
!
!         CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
!
          DO k = 1 , neleb
            iel = nelbor(k)
            xm( iel , cornsy(nulone(k),1) ) =
     &      xm( iel , cornsy(nulone(k),1) ) + xn(k)
            xm( iel , cornsy(nulone(k),2) ) =
     &      xm( iel , cornsy(nulone(k),2) ) + xn(k+nelebx)
          ENDDO
!
        ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
!
!         CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
!
          DO k = 1 , neleb
            iel = nelbor(k)
            xm( iel , cornsy(nulone(k),1) ) =
     &      xm( iel , cornsy(nulone(k),1) ) + xn(k)
            xm( iel , cornsy(nulone(k),2) ) =
     &      xm( iel , cornsy(nulone(k),2) ) + xn(k)
          ENDDO
!
        ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
!
!         CASE WHERE BOTH MATRICES ARE SYMMETRICAL
!
          DO k = 1 , neleb
            iel = nelbor(k)
            xm( iel , corsym(nulone(k)) ) =
     &      xm( iel , corsym(nulone(k)) ) + xn(k)
          ENDDO
!
        ELSE
          WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
99        FORMAT(1x,'OM1201 (BIEF) : TYPEXM = ',a1,' DOES NOT GO',
     &      /,1x,'FOR THE OPERATION : ',a8,' WITH TYPEXN = ',a1)
          CALL plante(1)
          stop
        ENDIF
!
!-----------------------------------------------------------------------
!
      ELSEIF(op(1:8).EQ.'M=M+TN  ') THEN
!
        CALL ovdb( 'X=X+Y   ' , dm , dn , z , c , nbor , nptfr )
!
        IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
!
!         CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
!
          DO k = 1 , neleb
            iel = nelbor(k)
            xm( iel , cornsy(nulone(k),1) ) =
     &      xm( iel , cornsy(nulone(k),1) ) + xn(k+nelebx)
            xm( iel , cornsy(nulone(k),2) ) =
     &      xm( iel , cornsy(nulone(k),2) ) + xn(k)
          ENDDO
!
        ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
!
!         CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
!
          DO k = 1 , neleb
            iel = nelbor(k)
            xm( iel , cornsy(nulone(k),1) ) =
     &      xm( iel , cornsy(nulone(k),1) ) + xn(k)
            xm( iel , cornsy(nulone(k),2) ) =
     &      xm( iel , cornsy(nulone(k),2) ) + xn(k)
          ENDDO
!
        ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
!
!         CASE WHERE BOTH MATRICES ARE SYMMETRICAL
!
          DO k = 1 , neleb
            iel = nelbor(k)
            xm( iel , corsym(nulone(k)) ) =
     &      xm( iel , corsym(nulone(k)) ) + xn(k)
          ENDDO
!
        ELSE
          WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
          CALL plante(1)
          stop
        ENDIF
!
!-----------------------------------------------------------------------
!
      ELSE
!
        WRITE(lu,141) op
141     FORMAT(1x,'OM1201 (BIEF) : UNKNOWN OPERATION : ',a8)
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END