om.f

Go to the documentation of this file.

!                   *************
                    SUBROUTINE om
!                   *************
!
     &( op , m , n , d , c , mesh )
!
!***********************************************************************
! BIEF   V7P0                                     21/08/2010
!***********************************************************************
!
!brief    OPERATIONS ON MATRICES.
!code
!+   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=N     '  : COPIES N IN M
!+      OP = 'M=CN    '  : MULTIPLIES N BY C
!+      OP = 'M=CM    '  : MULTIPLIES M BY C
!+      OP = 'M=M+CN  '  : ADDS CN TO M
!+      OP = 'M=MD    '  : M X D
!+      OP = 'M=DM    '  : D X M
!+      OP = 'M=DMD   '  : D X M X D
!+      OP = 'M=0     '  : SETS M TO 0 (TO BE CHECKED)
!+      OP = 'M=X(M)  '  : NOT SYMMETRICAL FORM OF M
!+      OP = 'M=TN    '  : COPIES TRANSPOSE OF N IN M
!+      OP = 'M=MSK(M)'  : MASKS M EXTRADIAGONAL TERMS
!
!warning  BE CAREFUL IF A NEW OPERATION IS ADDED TO THE LIST
!warning  IF OP CONTAINS N, IT THEN MEANS THAT MATRIX N IS USED
!
!history  ALGIANE FROEHLY
!+        13/02/2008
!+
!+   ADDED OM1113 AND OM1311
!
!history  J-M HERVOUET (LNHE)     ; F  LEPEINTRE (LNH)
!+        05/02/2010
!+        V6P0
!+   CALL TO OMSEGBOR MODIFIED, OMSEGPAR SUPPRESSED
!
!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  F. DECUNG (LNHE)
!+        01/01/2013
!+        V6P3
!+   omborseg added : operations on matrices with an edge-based storage
!         where N is a boundary matrix
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        13/03/2014
!+        V7P0
!+   Now written to enable different numbering of boundary points and
!+   boundary segments.
!history  R.NHEILI (Univerte de Perpignan, DALI)
!+        24/02/2016
!+        V7P3
!+      ADD MODASS=3
!
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| C              |-->| A GIVEN CONSTANT USED IN OPERATION OP
!| D              |-->| A DIAGONAL MATRIX
!| M              |<->| RESULTING MATRIX
!| MESH           |-->| MESH STRUCTURE
!| N              |-->| MATRIX USED IN FORMULA OP
!| OP             |-->| OPERATION TO BE DONE (SEE ABOVE)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_om => om
      USE declarations_telemac, ONLY : modass
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      CHARACTER(LEN=8), INTENT(IN)    :: OP
      TYPE(bief_obj)  , INTENT(INOUT), TARGET, OPTIONAL :: M
      TYPE(bief_obj)  , INTENT(IN), TARGET, OPTIONAL    :: N
      TYPE(bief_obj)  , INTENT(IN), TARGET, OPTIONAL    :: D
      DOUBLE PRECISION, INTENT(IN), OPTIONAL            :: C
      TYPE(bief_mesh) , INTENT(IN), OPTIONAL            :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELM1,IELM2,IELN1,IELN2,NELEM,NELMAX
      INTEGER NDIAGM,NDIAGN,NDIAGX,MSEG1,MSEG2,NSEG1,NSEG2
      INTEGER STOM,STON,NPTFR,MDIAGX
      INTEGER SIZXN,SZMXN,NETAGE
!
      LOGICAL YAN,YAD,YAC
      TYPE(bief_obj), POINTER :: NN
      TYPE(bief_obj), POINTER :: DD
      DOUBLE PRECISION CC
!
      CHARACTER(LEN=1) TYPDIM,TYPEXM,TYPDIN,TYPEXN
!
      INTEGER, DIMENSION(:), POINTER :: IKLE
!
      yan=.false.
      yad=.false.
      yac=.false.
      IF(inclus(op,'N')) yan=.true.
      IF(inclus(op,'D').OR.inclus(op,'MSK')) yad=.true.
      IF(inclus(op,'C')) yac=.true.
!
      IF(yan) THEN
        IF(PRESENT(n)) THEN
          nn=>n
        ELSE
          WRITE(lu,*) "---------------------------------"
          WRITE(lu,2) op
2         FORMAT(1x,'OM (BIEF) : N MISSING AND OPERATION ',a8,' ASKED')
          CALL plante(1)
          stop
        ENDIF
      ELSE
        nn=>m
      ENDIF
!
      IF(yad) THEN
        IF(PRESENT(d)) THEN
          dd=>d
        ELSE
          WRITE(lu,*) "---------------------------------"
          WRITE(lu,4) op
4         FORMAT(1x,'OM (BIEF) : D MISSING AND OPERATION ',a8,' ASKED')
          CALL plante(1)
          stop
        ENDIF
      ELSE
        dd=>m%D
      ENDIF
!
      IF(PRESENT(c)) THEN
        cc=c
      ELSE
        IF(yac) THEN
          WRITE(lu,*) "---------------------------------"
          WRITE(lu,6) op
6         FORMAT(1x,'OM (BIEF) : C MISSING AND OPERATION ',a8,' ASKED')
          CALL plante(1)
          stop
        ENDIF
      ENDIF
!
!-----------------------------------------------------------------------
!
!  CASE WHERE THE STRUCTURE OF M BECOMES THAT OF N
!
      IF(op(3:8).EQ.'N     '.OR.op(3:8).EQ.'CN    ') THEN
        CALL cpstmt(n,m)
      ELSEIF(op(3:8).EQ.'TN    ') THEN
        CALL cpstmt(n,m,trans=.true.)
      ENDIF
!
!  EXTRACTS THE CHARACTERISTICS OF MATRIX M
!
      typdim = m%TYPDIA
      typexm = m%TYPEXT
      stom = m%STO
      ndiagm = m%D%DIM1
      mdiagx = m%D%MAXDIM1
      ielm1 = m%ELMLIN
      ielm2 = m%ELMCOL
!
      IF(op(3:8).EQ.'X(M)  ') THEN
        IF(m%ELMLIN.NE.m%ELMCOL) THEN
          WRITE(lu,901) m%NAME
901       FORMAT(1x,'OM (BIEF) : M (REAL NAME: ',a6,') NOT SQUARE')
          WRITE(lu,701)
701       FORMAT(1x,'            IS ALREADY NON SYMMETRICAL')
          CALL plante(1)
          stop
        ENDIF
        IF(m%X%MAXDIM2*m%X%MAXDIM1.LT.
     &      bief_dim1_ext(ielm1,ielm2,stom,'Q',mesh)*
     &      bief_dim2_ext(ielm1,ielm2,stom,'Q',mesh)    ) THEN
            WRITE(lu,401) m%NAME
            WRITE(lu,801)
801         FORMAT(1x,'            TO BECOME NON SYMMETRICAL')
            CALL plante(1)
            stop
        ENDIF
      ENDIF
!
!-----------------------------------------------------------------------
!
!     EXTRACTS THE CHARACTERISTICS OF MATRIX N (OPTIONAL)
!
      IF(inclus(op,'N')) THEN
        typdin = nn%TYPDIA
        typexn = nn%TYPEXT
        ston   = nn%STO
        ndiagn = nn%D%DIM1
        ndiagx = nn%D%MAXDIM1
!
        sizxn = nn%X%DIM1
        szmxn = nn%X%MAXDIM1
!
! 07/02/03 : DIVISION BY NDIAGN=0 AVOIDED (SUBDOMAIN WITHOUT BOUNDARY POINTS
!            IN PARALLEL). COURTESY OLIVER GOETHEL (HANNOVER UNIVERSITY)
!
        IF(ndiagn.GT.0) THEN
          netage = ndiagm/ndiagn - 1
        ELSE
          netage = 0
        ENDIF
!
        ieln1 = nn%ELMLIN
        ieln2 = nn%ELMCOL
        IF(ndiagn.GT.mdiagx) THEN
          WRITE(lu,401) m%NAME
401       FORMAT(1x,'OM (BIEF) : M (REAL NAME: ',a6,') TOO SMALL')
          CALL plante(1)
          stop
        ENDIF
      ELSE
        ieln1 = ielm1
        ieln2 = ielm2
        ston = stom
      ENDIF
!
!-----------------------------------------------------------------------
!
!  DEPLOYMENT OF THE MESH STRUCTURE
!
!     STANDARD MATRIX
      IF(dimens(ielm1).EQ.mesh%DIM1) THEN
        ikle=>mesh%IKLE%I
        nelem = mesh%NELEM
        nelmax= mesh%NELMAX
      ELSE
!     BOUNDARY MATRIX
        ikle=>mesh%IKLBOR%I
        nelem  = mesh%NELEB
        nelmax = mesh%NELEBX
      ENDIF
!
      nptfr= mesh%NPTFR
!
!-----------------------------------------------------------------------
!
!  TRADITIONAL EBE STORAGE:
!
      IF(stom.EQ.1.AND.ston.EQ.1) THEN
!
      IF(ielm1.EQ.1.AND.ielm2.EQ.1) THEN
!
!     ELEMENTS WITH 2 POINTS
!
      CALL om0101(op,m%D%R,typdim,m%X%R,typexm,
     &            nn%D%R,typdin,nn%X%R,typexn,dd%R,cc,
     &            ikle,nelem,nelmax,ndiagm)
!
!
!     ELEMENTS WITH 3 POINTS
!
      ELSEIF( (ielm1.EQ.2 .AND.ielm2.EQ.2 ) .OR.
     &        (ielm1.EQ.11.AND.ielm2.EQ.11) .OR.
     &        (ielm1.EQ.61.AND.ielm2.EQ.61) .OR.
     &        (ielm1.EQ.81.AND.ielm2.EQ.81)       ) THEN
!
        IF( (ieln1.EQ.2 .AND.ieln2.EQ.2 ) .OR.
     &      (ieln1.EQ.11.AND.ieln2.EQ.11) .OR.
     &      (ieln1.EQ.61.AND.ieln2.EQ.61) .OR.
     &      (ieln1.EQ.81.AND.ieln2.EQ.81)         ) THEN
!
          IF (modass .EQ. 3) THEN
            CALL om1111(op , m%D%R,typdim,m%X%R,typexm ,
     &                  nn%D%R,typdin,nn%X%R,typexn,dd%R,cc,
     &                  ikle,nelem,nelmax,ndiagm,
     &                  dm_err=m%D%E, dn_err=nn%D%E,
     &                  d_err=d%E)
          ELSE
            CALL om1111(op , m%D%R,typdim,m%X%R,typexm ,
     &                  nn%D%R,typdin,nn%X%R,typexn,dd%R,cc,
     &                  ikle,nelem,nelmax,ndiagm)
          ENDIF
!
        ELSEIF(ieln1.EQ.1.AND.ieln2.EQ.1) THEN
!
          CALL om1101(op , m%D%R,typdim,m%X%R,typexm,
     &                     nn%D%R,typdin,nn%X%R,typexn,cc,
     &                     mesh%NULONE%I,mesh%NELBOR%I,
     &                     mesh%NBOR%I,
     &                     nelmax,ndiagm,ndiagn,mesh%NELEBX,
     &                     mesh%NELEB)
!
        ELSE
          WRITE(lu,101) m%NAME
          WRITE(lu,201) ielm1,ielm2
          WRITE(lu,151) nn%NAME
          WRITE(lu,251) ieln1,ieln2
          WRITE(lu,301)
          CALL plante(1)
          stop
        ENDIF
!
!     ELEMENTS WITH 4 POINTS
!
      ELSEIF( (ielm1.EQ.21.AND.ielm2.EQ.21) .OR.
     &        (ielm1.EQ.71.AND.ielm2.EQ.71) .OR.
     &        (ielm1.EQ.31.AND.ielm2.EQ.31) .OR.
     &        (ielm1.EQ.51.AND.ielm2.EQ.51) .OR.
     &        (ielm1.EQ.12.AND.ielm2.EQ.12)      ) THEN
!
        IF(   (ieln1.EQ.21.AND.ieln2.EQ.21) .OR.
     &        (ieln1.EQ.71.AND.ieln2.EQ.71) .OR.
     &        (ieln1.EQ.31.AND.ieln2.EQ.31) .OR.
     &        (ieln1.EQ.51.AND.ieln2.EQ.51) .OR.
     &        (ieln1.EQ.12.AND.ieln2.EQ.12)      ) THEN
!
          CALL om2121(op , m%D%R,typdim,m%X%R,typexm ,
     &                nn%D%R,typdin,nn%X%R,typexn,dd%R,cc,
     &                ikle,nelem,nelmax,ndiagm)
        ELSEIF(  (ielm1.EQ.12.AND.ielm2.EQ.12) .AND.
     &           (ieln1.EQ.1 .AND.ieln2.EQ.1 )   ) THEN
!
          CALL om1201(op , m%D%R,typdim,m%X%R,typexm,
     &                     nn%D%R,typdin,nn%X%R,typexn,cc,
     &                     mesh%NULONE%I,mesh%NELBOR%I,
     &                     mesh%NBOR%I,
     &                     nelmax,ndiagm,ndiagn,mesh%NELEBX,
     &                     mesh%NELEB)
        ELSEIF(( (ielm1.EQ.51.AND.ielm2.EQ.51) .AND.
     &           (ieln1.EQ.61.AND.ieln2.EQ.61))   ) THEN
!         PRISMS SPLIT IN TETRAHEDRONS M INTERIOR MATRIX
!                                      N LATERAL BOUNDARY MATRIX
          CALL om5161(op,m%D%R,typdim,m%X%R,typexm,
     &                nn%D%R,typdin,nn%X%R,typexn,cc,
     &                mesh%NULONE%I,mesh%NELBOR%I,mesh%NBOR%I,
     &                nelmax,ndiagn,sizxn,szmxn)
!
        ELSEIF(( (ielm1.EQ.31.AND.ielm2.EQ.31) .AND.
     &           (ieln1.EQ.81.AND.ieln2.EQ.81))   ) THEN
!         NOT STRUCTURED TETRAHEDRONS    M INTERIOR MATRIX
!                                        N LATERAL BOUNDARY MATRIX
          CALL om3181(op,m%D%R,typdim,m%X%R,typexm,
     &                nn%D%R,typdin,nn%X%R,typexn,cc,
     &                mesh%NULONE%I,mesh%NELBOR%I,mesh%NBOR%I,
     &                nelmax,ndiagn,mesh%NELEB)
!
        ELSEIF(  (ielm1.EQ.51.AND.ielm2.EQ.51) .AND.
     &           (ieln1.EQ.11.AND.ieln2.EQ.11)   ) THEN
!         PRISMS SPLIT IN TETRAHEDRONS M INTERIOR MATRIX
!                                      N BOTTOM (NB) OR SURFACE (NS) BOUNDARY MATRIX
!                                      OPERATIONS M+NB AND M+NS
          CALL om5111(op , m%D%R,typdim,m%X%R,typexm,
     &                nn%D%R,typdin,nn%X%R,typexn,
     &                ndiagn,ndiagx,sizxn,netage,nelmax)
!
        ELSE
          WRITE(lu,101) m%NAME
          WRITE(lu,201) ielm1,ielm2
          WRITE(lu,151) nn%NAME
          WRITE(lu,251) ieln1,ieln2
          WRITE(lu,301)
          CALL plante(1)
          stop
        ENDIF
!
!     RECTANGULAR MATRICES WITH 3 AND 4 POINTS
!
      ELSEIF(ielm1.EQ.11.AND.ielm2.EQ.12) THEN
!
        IF((ieln1.EQ.11.AND.ieln2.EQ.12).OR.
     &     (ieln1.EQ.12.AND.ieln2.EQ.11.AND.op(1:4).EQ.'M=TN')) THEN
          CALL om1112(op , m%D%R,typdim,m%X%R,typexm ,
     &                     nn%D%R,typdin,nn%X%R,typexn,dd%R,cc,
     &                     ikle,nelem,nelmax,ndiagm)
        ELSE
          WRITE(lu,101) m%NAME
          WRITE(lu,201) ielm1,ielm2
          WRITE(lu,151) nn%NAME
          WRITE(lu,251) ieln1,ieln2
          WRITE(lu,301)
          CALL plante(1)
          stop
        ENDIF
!
!     RECTANGULAR MATRICES WITH 4 AND 3 POINTS
!
      ELSEIF(ielm1.EQ.12.AND.ielm2.EQ.11) THEN
!
        IF((ieln1.EQ.12.AND.ieln2.EQ.11).OR.
     &     (ieln1.EQ.11.AND.ieln2.EQ.12.AND.op(1:4).EQ.'M=TN')) THEN
          CALL om1211(op , m%D%R,typdim,m%X%R,typexm ,
     &                nn%D%R,typdin,nn%X%R,typexn,dd%R,cc,
     &                ikle,nelem,nelmax,ndiagm)
        ELSE
          WRITE(lu,101) m%NAME
          WRITE(lu,201) ielm1,ielm2
          WRITE(lu,151) nn%NAME
          WRITE(lu,251) ieln1,ieln2
          WRITE(lu,301)
          CALL plante(1)
          stop
        ENDIF
!
!     ELEMENTS WITH 6 POINTS
!
      ELSEIF( (ielm1.EQ.41.AND.ielm2.EQ.41).OR.
     &        (ielm1.EQ.13.AND.ielm2.EQ.13)      ) THEN
!
        IF( (ieln1.EQ.41.AND.ieln2.EQ.41).OR.
     &      (ieln1.EQ.13.AND.ieln2.EQ.13)        ) THEN
!
          CALL om4141(op,m%D%R,typdim,m%X%R,typexm,
     &                nn%D%R,typdin,nn%X%R,typexn,dd%R,cc,
     &                ikle,nelem,nelmax,ndiagm,m%STOX)
!
        ELSEIF(ieln1.EQ.71.AND.ieln2.EQ.71) THEN
!
          CALL om4121(op,m%D%R,typdim,m%X%R,typexm,
     &                   nn%D%R,typdin,nn%X%R,typexn,cc,
     &                   mesh%NULONE%I,mesh%NELBOR%I,mesh%NBOR%I,
     &                   nelmax,ndiagn,sizxn,szmxn)
!
        ELSEIF(ieln1.EQ.11.AND.ieln2.EQ.11) THEN
!
          CALL om4111(op , m%D%R,typdim,m%X%R,typexm,
     &                     nn%D%R,typdin,nn%X%R,typexn,
     &                     ndiagn,ndiagx,sizxn,netage,nelmax)
!
        ELSEIF(  (ielm1.EQ.13.AND.ielm2.EQ.13) .AND.
     &           (ieln1.EQ.2 .AND.ieln2.EQ.2 )   ) THEN
!
          CALL om1302(op , m%D%R,typdim,m%X%R,typexm,
     &                     nn%D%R,typdin,nn%X%R,typexn,cc,
     &                     mesh%NULONE%I,mesh%NELBOR%I,
     &                     mesh%NBOR%I,
     &                     nelmax,ndiagm,nptfr,mesh%NELEBX,
     &                     mesh%NELEB)
!
        ELSE
          WRITE(lu,101) m%NAME
          WRITE(lu,201) ielm1,ielm2
          WRITE(lu,151) nn%NAME
          WRITE(lu,251) ieln1,ieln2
          WRITE(lu,301)
          CALL plante(1)
          stop
        ENDIF
!
!     RECTANGULAR MATRICES WITH 3 AND 6 POINTS
!
      ELSEIF(ielm1.EQ.11.AND.ielm2.EQ.13) THEN
!
        IF((ieln1.EQ.11.AND.ieln2.EQ.13).OR.
     &     (ieln1.EQ.13.AND.ieln2.EQ.11.AND.op(1:4).EQ.'M=TN')) THEN
          CALL om1113(op , m%D%R,typdim,m%X%R,typexm ,
     &                     nn%D%R,typdin,nn%X%R,typexn,dd%R,cc,
     &                     ikle,nelem,nelmax,ndiagm)
        ELSE
          WRITE(lu,101) m%NAME
          WRITE(lu,201) ielm1,ielm2
          WRITE(lu,151) nn%NAME
          WRITE(lu,251) ieln1,ieln2
          WRITE(lu,301)
          CALL plante(1)
          stop
        ENDIF
!
!     RECTANGULAR MATRICES WITH 6 AND 3 POINTS
!
      ELSEIF(ielm1.EQ.13.AND.ielm2.EQ.11) THEN
!
        IF((ieln1.EQ.13.AND.ieln2.EQ.11).OR.
     &     (ieln1.EQ.11.AND.ieln2.EQ.13.AND.op(1:4).EQ.'M=TN')) THEN
          CALL om1311(op , m%D%R,typdim,m%X%R,typexm ,
     &                nn%D%R,typdin,nn%X%R,typexn,dd%R,cc,
     &                ikle,nelem,nelmax,ndiagm)
        ELSE
          WRITE(lu,101) m%NAME
          WRITE(lu,201) ielm1,ielm2
          WRITE(lu,151) nn%NAME
          WRITE(lu,251) ieln1,ieln2
          WRITE(lu,301)
          CALL plante(1)
          stop
        ENDIF
!
!  COMBINATION OF IELM1 AND IELM2 NOT IMPLEMENTED: ERROR
!
      ELSE
        WRITE(lu,101) m%NAME
        WRITE(lu,201) ielm1,ielm2
        WRITE(lu,411) stom,ston
        WRITE(lu,301)
        CALL plante(1)
        stop
      ENDIF
!
      ELSEIF(stom.EQ.3.AND.ston.EQ.3) THEN
!
!  STORAGE BY SEGMENT
!
        IF(m%ELMCOL.NE.nn%ELMCOL.OR.m%ELMLIN.NE.nn%ELMLIN) THEN
!         WRITE(LU,*) 'M ET N DE STRUCTURES DIFFERENTES',M%ELMCOL
!         CALL PLANTE(1)
!         STOP
!
!     EDGE-BASED STORAGE FOR M AND N
!     THIS CAN HAPPEN ONLY WHEN N IS A BOUNDARY MATRIX (FD : REALLY?)
!     TESTED SO FAR FOR TETRA (81) AND TRIANGLES (31)
!     BORDER SEGMENTS ARE LINKED TO THE NSEGBOR FIRST SEGMENTS
!
          mseg1 = bief_nbseg(m%ELMLIN,mesh)
          mseg2 = bief_nbseg(m%ELMCOL,mesh)
          nseg1 = bief_nbseg(nn%ELMLIN,mesh)
          nseg2 = bief_nbseg(nn%ELMCOL,mesh)
!
          CALL omborseg(op,m%D%R,m%X%R,typexm,
     &         nn%D%R,nn%X%R,typexn,cc,
     &         ndiagn,mseg1,mseg2,nseg1,nseg2,
     &         mesh%NBOR%I)
!
        ELSE
!
          nseg1 = bief_nbseg(m%ELMLIN,mesh)
          nseg2 = bief_nbseg(m%ELMCOL,mesh)
!
!         IN LINEAR-QUADRATIC RECTANGULAR MATRICES, PURELY QUADRATIC
!         SEGMENTS ARE NOT CONSIDERED (NUMBER 13,14 AND 15, SO 3 PER ELEMENT)
!
          IF(ielm1.EQ.11.AND.ielm2.EQ.13) THEN
            nseg2=nseg2-3*nelem
          ELSEIF(ielm1.EQ.13.AND.ielm2.EQ.11) THEN
            nseg1=nseg1-3*nelem
          ENDIF
!
!         IN LINEAR-QUADRATIC RECTANGULAR MATRICES, PURELY QUADRATIC
!         SEGMENTS ARE NOT CONSIDERED (NUMBER 13,14 AND 15, SO 3 PER ELEMENT)
!
          CALL omseg(op , m%D%R,typdim,m%X%R,typexm,
     &                    nn%D%R,typdin,nn%X%R,typexn,dd%R,cc,
     &                    ndiagm,nseg1,nseg2,mesh%GLOSEG%I,
     &                    mesh%GLOSEG%MAXDIM1)
!
        ENDIF
!
      ELSEIF(stom.EQ.3.AND.ston.EQ.1) THEN
!
!       EDGE-BASED STORAGE FOR M AND EBE FOR N
!       THIS CAN HAPPEN ONLY WHEN N IS A BOUNDARY MATRIX
!
        IF(  ( m%ELMLIN.EQ.11.AND. m%ELMCOL.EQ.11.AND.
     &        nn%ELMLIN.EQ.1 .AND.nn%ELMCOL.EQ.1) .OR.
     &       ( m%ELMLIN.EQ.12.AND. m%ELMCOL.EQ.12.AND.
     &        nn%ELMLIN.EQ.1 .AND.nn%ELMCOL.EQ.1) .OR.
     &       ( m%ELMLIN.EQ.13.AND. m%ELMCOL.EQ.13.AND.
     &        nn%ELMLIN.EQ.1 .AND.nn%ELMCOL.EQ.1) .OR.
     &       ( m%ELMLIN.EQ.13.AND. m%ELMCOL.EQ.13.AND.
     &        nn%ELMLIN.EQ.2 .AND.nn%ELMCOL.EQ.2)      ) THEN
!
          nseg1 = bief_nbseg(m%ELMLIN,mesh)
          nseg2 = bief_nbseg(m%ELMCOL,mesh)
!
          CALL omsegbor(op,m%D%R,typdim,m%X%R,typexm,
     &                     nn%D%R,typdin,nn%X%R,typexn,cc,
     &                     ndiagm,nseg1,mesh%NBOR%I,
     &                     nptfr,m%ELMLIN,nn%ELMLIN,
     &                     bief_nbseg(11,mesh),
     &                     mesh%IKLBOR%I,mesh%NELEBX,mesh%NELEB)
!
        ELSE
          WRITE(lu,*) 'OM : UNEXPECTED CASE IN SEGMENT STORAGE'
          WRITE(lu,*) '     M%ELMLIN=',m%ELMLIN
          WRITE(lu,*) '     M%ELMCOL=',m%ELMCOL
          WRITE(lu,*) '     M%NAME=',m%NAME
          WRITE(lu,*) '     NN%ELMLIN=',nn%ELMLIN
          WRITE(lu,*) '     NN%ELMCOL=',nn%ELMCOL
          WRITE(lu,*) '     NN%NAME=',nn%NAME
          WRITE(lu,*) '     IMPLEMENTATION MISSING'
          CALL plante(1)
          stop
        ENDIF
!
!  COMBINATION OF IELM1 AND IELM2 NOT IMPLEMENTED: ERROR
!
!     ELSE
!        WRITE(LU,101) M%NAME
!        WRITE(LU,201) IELM1,IELM2
!        WRITE(LU,411) STOM,STON
!        WRITE(LU,301)
!        CALL PLANTE(1)
!        STOP
!     ENDIF
!
      ELSE
!
!  STORAGE COMBINATION NOT IMPLEMENTED
!
        WRITE(lu,101) m%NAME
        WRITE(lu,411) stom,ston
        WRITE(lu,301)
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!  RE-ENCODES THE NEW TYPE
!
      m%TYPDIA = typdim
      m%TYPEXT = typexm
      IF(op(3:8).EQ.'X(M)  ') THEN
        m%X%DIM1=bief_dim1_ext(ielm1,ielm2,stom,'Q',mesh)
        m%X%DIM2=bief_dim2_ext(ielm1,ielm2,stom,'Q',mesh)
      ENDIF
!
!-----------------------------------------------------------------------
!
101   FORMAT(1x,'OM (BIEF) : MATRIX  M (REAL NAME:',a6,')')
151   FORMAT(1x,'OM (BIEF) : MATRIX  N (REAL NAME:',a6,')')
201   FORMAT(1x,'            IELM1 = ',1i6,' IELM2 = ',1i6)
251   FORMAT(1x,'            IELN1 = ',1i6,' IELN2 = ',1i6)
301   FORMAT(1x,'            THIS CASE IS NOT IMPLEMENTED')
411   FORMAT(1x,'AND STORAGES   M  : ',1i6,' STON  : ',1i6)
!
!-----------------------------------------------------------------------
!
      RETURN
      END