matrix.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE matrix
!                   *****************
!
     &(m,op,formul,ielm1,ielm2,xmul,f,g,h,u,v,w,mesh,msk,maskel)
!
!***********************************************************************
! BIEF   V6P3                                   21/08/2010
!***********************************************************************
!
!brief    OPERATIONS BETWEEN MATRICES.
!+
!+            THE MATRIX IS IDENTIFIED BY THE FORMULATION IN
!+                CHARACTER STRING FORMUL.
!code
!+     OP IS A STRING OF 8 CHARACTERS, WHICH INDICATES HOW M IS
!+     MODIFIED. THE SYNTAX IS THE SAME AS THAT OF OM, FOR EXAMPLE.
!+
!+     OP='M=N     '
!+     OP='M=TN    '
!+     OP='M=M+N   '
!+     OP='M=M+TN  '
!+
!+     ALL THE OPERATIONS IN OM WHICH HAVE N ON THE RIGHT OF THE
!+     = SIGN ARE VALID.
!
!code
!+-----------------------------------------------------------------------
!+  MEANING OF IELM AND IELM2
!+
!+  TYPE OF ELEMENT      NUMBER OF POINTS
!+
!+  10 : P0 TRIANGLE            1
!+  11 : P1 TRIANGLE            3
!+  12 : QUASI-BUBBLE TRIANGLE  4
!+
!+  20 : Q0 QUADRILATERAL       1
!+  21 : Q1 QUADRILATERAL       4
!+
!+  40 : TELEMAC-3D P0 PRISMS   1
!+  41 : TELEMAC-3D P1 PRISMS   6
!+
!+-----------------------------------------------------------------------
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        25/06/2008
!+
!+   DOES NOT CALL MATRIY IF NUMBER OF ELEMENTS IS 0
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        14/10/2009
!+        V6P0
!+   ARGUMENTS ADDED TO ASSEX3
!
!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 R&D, LNHE)
!+        11/01/2013
!+        V6P3
!+   Arguments added to MATRIY
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        13/03/2014
!+        V7P0
!+   Cases where there is no boundary element secured.
!
!history  R.NHEILI (Univerte de Perpignan, DALI)
!+        24/02/2016
!+        V7P3
!+      ADD MODASS=3
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| F              |-->| FUNCTION USED IN THE FORMULA
!| FORMUL         |-->| FORMULA DESCRIBING THE RESULTING MATRIX
!| G              |-->| FUNCTION USED IN THE FORMULA
!| H              |-->| FUNCTION USED IN THE FORMULA
!| IELM1          |-->| TYPE OF ELEMENT FOR LINES
!| IELM2          |-->| TYPE OF ELEMENT FOR COLUMNS
!| M              |<->| MATRIX TO BE BUILT OR MODIFIED
!| MASKEL         |-->| MASKING OF ELEMENTS
!|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
!| MESH           |-->| MESH STRUCTURE
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
!| OP             |-->| OPERATION TO BE DONE. IF '=' THEN M=...
!|                |   |                       IF '+' THEN M=M+
!| U              |-->| FUNCTION USED IN THE FORMULA (COMPONENT OF VECTOR)
!| V              |-->| FUNCTION USED IN THE FORMULA (COMPONENT OF VECTOR)
!| W              |-->| FUNCTION USED IN THE FORMULA (COMPONENT OF VECTOR)
!| XMUL           |-->| COEFFICIENT FOR MULTIPLICATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_matrix => matrix
      USE declarations_telemac, ONLY : modass
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)            :: IELM1,IELM2
!
      DOUBLE PRECISION, INTENT(IN)   :: XMUL
!
      LOGICAL, INTENT(IN)            :: MSK
!
      CHARACTER(LEN=16), INTENT(IN)  :: FORMUL
      CHARACTER(LEN=8), INTENT(IN)   :: OP
!
      TYPE(bief_obj), INTENT(IN)     :: F,G,H,U,V,W,MASKEL
      TYPE(bief_obj), INTENT(INOUT)  :: M
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER NELMAX,NELEM,NPT,SS,I
      INTEGER, DIMENSION(:), POINTER :: IKLE
      INTEGER, DIMENSION(:), POINTER :: ELTSEG,ORISEG
      DOUBLE PRECISION C
      LOGICAL LEGO
!
!-----------------------------------------------------------------------
!
!     STORES 1 FOR THE WROKING ARRAY
!     CAN BE MODIFIED BY ASSEXT THEREAFTER
!
      mesh%M%STO = 1
!
!-----------------------------------------------------------------------
!  CALLS THE SHUNTING AND ASSEMBLY SUBROUTINE : MATRIY
!-----------------------------------------------------------------------
!
!     LEGO CAN BE MODIFIED BY MATRIY
      lego = .true.
!
      IF(dimens(ielm1).EQ.mesh%DIM1) THEN
!       NORMAL MATRIX
        nelem  = mesh%NELEM
        nelmax = mesh%NELMAX
        ikle   =>mesh%IKLE%I
        eltseg=>mesh%ELTSEG%I
        oriseg=>mesh%ORISEG%I
        IF (modass .EQ.3) THEN
          DO i=1, mesh%M%D%DIM1
            mesh%M%D%E=0.d0
          ENDDO
        ENDIF
      ELSE
!       BOUNDARY MATRIX
        nelem  = mesh%NELEB
        nelmax = mesh%NELEBX
        ikle   =>mesh%IKLBOR%I
        eltseg=>mesh%ELTSEGBOR%I
        oriseg=>mesh%ORISEGBOR%I
      ENDIF
!
!     MATRIY FILLS THE DIAGONAL AND EXTRA DIAGONAL TERMS
!
!     REFLECTS CHOICE OF PRE-ASSEMBLY STORAGE
      IF(m%STO.EQ.1.OR.m%STO.EQ.3) THEN
        ss = 1
      ELSEIF(m%STO.EQ.2) THEN
        ss = 2
      ELSE
        WRITE(lu,*) 'UNKNOWN STORAGE IN MATRIX'
        CALL plante(1)
        stop
      ENDIF
!
      IF(nelem.GT.0) THEN
        CALL matriy(formul,mesh%M%X%R,
     &              mesh%M%TYPDIA,mesh%M%TYPEXT,
     &              xmul,f,g,h,u,v,w,
     &              f%R,g%R,h%R,u%R,v%R,w%R,
     &              mesh%W%R,lego,
     &              mesh%XEL%R,mesh%YEL%R,mesh%ZEL%R,
     &              mesh%X%R,mesh%Y%R,mesh%Z%R,
     &              mesh%SURFAC%R,mesh%LGSEG%R,
     &              mesh%IKLE%I,mesh%IKLBOR%I,mesh%NBOR%I,
     &              mesh%NELBOR%I,mesh%NULONE%I,
     &              mesh%NELEM,mesh%NELMAX,
     &              mesh%NELEB,mesh%NELEBX,ielm1,ielm2,ss,
     &              mesh%NPOIN/bief_nbpts(11,mesh),mesh,nelmax,
     &              mesh%M%STOX)
      ENDIF
!
!  UPDATES THE INFORMATION OF THE MATRIX
!
      npt = bief_nbpts(min(ielm1,ielm2),mesh)
!
      IF(npt.GT.mesh%M%D%MAXDIM1) THEN
        WRITE(lu,501) mesh%M%NAME
        WRITE(lu,2001) ielm1
        WRITE(lu,3001) ielm2
        CALL plante(1)
        stop
      ENDIF
!
      mesh%M%D%ELM  = min(ielm1,ielm2)
      mesh%M%D%DIM1 = npt
      mesh%M%ELMLIN = ielm1
      mesh%M%ELMCOL = ielm2
!
!  ASSEMBLES THE DIAGONAL (POSSIBLY)
!
      IF(lego.AND.mesh%M%TYPDIA.EQ.'Q'.AND.nelem.GT.0) THEN
!
        IF (modass .EQ.3) THEN
          CALL assvec(mesh%M%D%R,
     &                  ikle,npt,nelem,nelmax,
     &                  mesh%W%R,lego,mesh%LV,msk,maskel%R,
     &                  bief_nbpel(mesh%M%D%ELM,mesh),mesh%M%D%E)
!
        ELSE
          CALL assvec(mesh%M%D%R,
     &                  ikle,npt,nelem,nelmax,
     &                  mesh%W%R,lego,mesh%LV,msk,maskel%R,
     &                  bief_nbpel(mesh%M%D%ELM,mesh))
        ENDIF
!
      ENDIF
!
!  MASKS EXTRA-DIAGONAL TERMS (POSSIBLY)
!
      IF(msk.AND.nelem.GT.0) THEN
        CALL om('M=MSK(M)',m=mesh%M,n=mesh%M,d=maskel,mesh=mesh)
      ENDIF
!
!  ASSEMBLES EXTRA-DIAGONAL TERMS (POSSIBLY)
!
      IF(m%STO.EQ.3) THEN
!       COPIES THE DIAGONAL
        CALL os('X=Y     ',x=mesh%MSEG%D,y=mesh%M%D)
        mesh%MSEG%TYPDIA(1:1)='Q'
!       ASSEMBLES EXTRA-DIAGONAL TERMS
        IF(mesh%M%TYPEXT.EQ.'Q'.OR.mesh%M%TYPEXT.EQ.'S') THEN
!
!         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!         CASE OF MATRICES WITH INVERTED STORAGE OF OFF-DIAGONAL TERMS
!         SEE INVERSION OF DIM1_EXT AND DIM2_EXT IN CALL TO ASSEX3
!
!         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
          IF(mesh%M%STOX.EQ.2) THEN
            CALL assex3(mesh%MSEG%X%R,mesh%M%STO,mesh%M%NAME,
     &                  mesh%M%ELMLIN,mesh%M%ELMCOL,
     &                  mesh%M%TYPEXT,mesh%M%X%R,
     &                  bief_dim2_ext(ielm1,ielm2,mesh%M%STO,
     &                                mesh%M%TYPEXT,mesh),
     &                  bief_dim1_ext(ielm1,ielm2,mesh%M%STO,
     &                                mesh%M%TYPEXT,mesh),
     &                  mesh%M%STOX,
     &                  mesh,nelmax,eltseg,oriseg)
          ELSEIF(mesh%M%STOX.EQ.1) THEN
            CALL assex3(mesh%MSEG%X%R,mesh%M%STO,mesh%M%NAME,
     &                  mesh%M%ELMLIN,mesh%M%ELMCOL,
     &                  mesh%M%TYPEXT,mesh%M%X%R,
     &                  bief_dim1_ext(ielm1,ielm2,mesh%M%STO,
     &                                mesh%M%TYPEXT,mesh),
     &                  bief_dim2_ext(ielm1,ielm2,mesh%M%STO,
     &                                mesh%M%TYPEXT,mesh),
     &                  mesh%M%STOX,
     &                  mesh,nelmax,eltseg,oriseg)
          ELSE
            WRITE(lu,*) 'MATRIX: UNKNOWN STORAGE:'
            WRITE(lu,*) 'MESH%M%STOX= : ',mesh%M%STOX
            CALL plante(1)
            stop
          ENDIF
        ENDIF
        mesh%MSEG%TYPEXT=mesh%M%TYPEXT
        mesh%MSEG%ELMLIN = ielm1
        mesh%MSEG%ELMCOL = ielm2
        mesh%MSEG%D%ELM  = min(ielm1,ielm2)
        mesh%MSEG%D%DIM1 = npt
        mesh%MSEG%X%DIM1 = bief_dim1_ext(ielm1,ielm2,m%STO,
     &                                   mesh%MSEG%TYPEXT,mesh)
        mesh%MSEG%X%DIM2 = bief_dim2_ext(ielm1,ielm2,m%STO,
     &                                   mesh%MSEG%TYPEXT,mesh)
      ENDIF
!
!     DIMENSIONS OF THE ARRAY WITH EXTRADIAGONAL TERMS
!     BEWARE M%STO (NOT MESH%M%STO BECAUSE IT EQUALS 1)
!                   SEE BEGINNING OF SUBROUTINE
!
      mesh%M%X%DIM1 = bief_dim1_ext(ielm1,ielm2,m%STO,
     &                              mesh%M%TYPEXT,mesh)
      mesh%M%X%DIM2 = bief_dim2_ext(ielm1,ielm2,m%STO,
     &                              mesh%M%TYPEXT,mesh)
!
!-----------------------------------------------------------------------
!  UPDATES M AFTER WORK ON MESH%M IS COMPLETE
!-----------------------------------------------------------------------
!
      IF(nelem.GT.0) THEN
!
        IF(m%STO.EQ.1) THEN
          CALL om(op, m, mesh%M, f, c, mesh)
        ELSEIF(m%STO.EQ.3) THEN
          CALL om(op, m, mesh%MSEG, f, c, mesh)
        ELSE
          WRITE(lu,*) 'MATRIX: UNKNOWN STORAGE : ',m%STO
          CALL plante(1)
          stop
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
!
501   FORMAT(1x,'MATRIX (BIEF) : MATRIX ',a6,' TOO SMALL')
2001  FORMAT(1x,'                FOR IELM1 = ',1i6)
3001  FORMAT(1x,'                AND IELM2 = ',1i6)
!
!-----------------------------------------------------------------------
!
      RETURN
      END