centre_mass_seg.f

Go to the documentation of this file.

!                       **************************
                        SUBROUTINE centre_mass_seg
!                       **************************
!
     &(x,y,coord_g,ikle,npoin,eltseg,oriseg,nelem,nseg,jmi,cmi,gloseg,
     & ifabor,mesh)
!
!***********************************************************************
! BIEF   V6P3                                                18/12/2012
!***********************************************************************
!
!brief    GIVES COORDINATES OF CENTRE OF GRAVITY OF TRIANGLES RIGHT AND
!           LEFT OF SEGMENTS I.E. FOR A SEGMENT ISEG:
!            - COORD_G(1,ISEG) AND COORD_G(2,ISEG) ARE THE COORDINATES OF
!              THE CENTRE OF GRAVITY OF THE FIRST NEIGHBORING TRIANGLE
!            - COORD_G(3,ISEG) AND COORD_G(4,ISEG) ARE THE COORDINATES OF
!              THE CENTRE OF GRAVITY OF THE SECOND NEIGHBORING TRIANGLE
!         GIVES COORDINATES OF THE MIDDLE POINT OF SEGMENT G1G2 (STOCKED
!             AT CMI(1,ISEG) AND CMI(2,ISEG)) AND THE ELEMENT NUMBER TO
!             TO WHICH BELONGS CMI
!
!history  R.ATA
!+        18/12/12
!+        V6P3
!+
!history  R. ATA & J-M HERVOUET
!+        19/11/2013
!+        V6P3
!+    Optimisation and simplification.
!
!history  R. ATA & J-M HERVOUET
!+        01/12/2013
!+        V6P3
!+    add verification tests
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CMI            |<--| COORDINATES OF MID-INTERFACE POINTS
!| COORD_G        |<--| CENTER OF MASS OF ELEMENTS NEIGHBORS OF AN EDGE
!|                |   |  COORD_G(1,ISEG) AND COORD_G(2,ISEG) ARE X1 AND Y1
!|                |   |  COORD_G(3,ISEG) AND COORD_G(4,ISEG) ARE X2 AND Y2
!| ELTSEG         |-->| SEGMENTS FORMING AN ELEMENT
!| GLOSEG         |-->|  GLOBAL NUMBERS OF VERTICES OF SEGMENTS
!| IKLE           |-->| CONNECTIVITY TABLE.
!| IFABOR         |-->| IFABOR(IEL,I) IS THE ELEMENT BEHIND THE EDGE I OF
!|                |---|  ELEMENT IEL, OTHERWISE
!|                |---|     IFABOR(IEL,I) = -2 : THIS IS INTERFACE EDGE
!|                |---|     IFABOR(IEL,I) = 0  : THIS IS BOUNDARY EDGE
!|                |---|     IFABOR(IEL,I) = -1 : THIS IS LIQUID BOUNDARY EDGE
!| JMI            |<--| NUMBER OF TRIANGLE TO WHICH BELONGS THE
!|                |   |       MID-INTERFACE POINT.
!| MESH           |-->| MESH STRUCTURE
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| NSEG           |-->| NUMBER OF SEGMENTS
!|                |---|
!| ORISEG         |-->| ORIENTATION OF SEGMENTS FORMING AN
!|                |   |        ELEMENT 1:ANTI 2:CLOCKWISE
!| X              |-->| ABSCISSAE OF POINTS IN THE MESH
!| Y              |-->| ORDINATES OF POINTS IN THE MESH
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ONLY : inpoly
      USE bief_def
      IMPLICIT  NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NSEG,NELEM,NPOIN
      INTEGER, INTENT(IN)             :: IKLE(nelem,3)
      INTEGER, INTENT(IN)             :: ELTSEG(nelem,3)
      INTEGER, INTENT(IN)             :: ORISEG(nelem,3)
      INTEGER, INTENT(INOUT)          :: JMI(nseg)
      INTEGER, INTENT(IN)             :: GLOSEG(nseg,2)
      DOUBLE PRECISION, INTENT(INOUT) :: CMI(2,nseg)
      DOUBLE PRECISION, INTENT(INOUT) :: COORD_G(nseg,4)
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin),Y(npoin)
      INTEGER, INTENT(IN)             :: IFABOR(nelem,3)
      TYPE(bief_mesh), INTENT(INOUT)  :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER ISEG,NB1,NB2,IELEM,I,I1,I2,I3
      DOUBLE PRECISION XEL,YEL,XG1,XG2,YG1,YG2,XSOM(3),YSOM(3)
      DOUBLE PRECISION X_MIDPOINT,Y_MIDPOINT,DET
      DOUBLE PRECISION, ALLOCATABLE :: TMP_CMI1(:), TMP_CMI2(:)
!
!-----------------------------------------------------------------------
!   INITIALIZATION OF VARIABLES
!
      DO i=1,nseg
        coord_g(i,1) = 0.d0
        coord_g(i,2) = 0.d0
        coord_g(i,3) = 0.d0
        coord_g(i,4) = 0.d0
        cmi(1,i)     = 0.d0
        cmi(2,i)     = 0.d0
        jmi(i)       = 0
      ENDDO
      IF(ncsize.GT.1)THEN
        DO i=1,nseg
          mesh%MSEG%X%R(i)=0.d0
        ENDDO
      ENDIF
!    RECUPERATE THE COORDINATES OF CENTER OF GRAVITY OF THE TRIANGLES
!    AT THE RIGHT AND AT THE LEFT OF THE SEGMENT BETWEEN TWO NODES
!
      DO ielem=1, nelem
        i1 = ikle(ielem,1)
        i2 = ikle(ielem,2)
        i3 = ikle(ielem,3)
        xel = (x(i1)+x(i2)+x(i3))/3.0d0
        yel = (y(i1)+y(i2)+y(i3))/3.0d0
        DO i = 1,3
          iseg = eltseg(ielem,i)
          IF(oriseg(ielem,i).EQ.1)THEN ! G IS LEFT OF THE EDGE
            coord_g(iseg,1) =xel
            coord_g(iseg,2) =yel
          ELSEIF(oriseg(ielem,i).EQ.2)THEN !G IS RIGHT OF THE EDGE
            coord_g(iseg,3) =xel
            coord_g(iseg,4) =yel
          ENDIF
        ENDDO
      ENDDO
!
      IF(ncsize.GT.1) THEN
        CALL parcom2_seg(coord_g(1,1),
     &                   coord_g(1,2),
     &                   coord_g(1,2), ! NO EFFECT FOR THIS ONE
     &                   nseg,1,1,2,mesh,1,11)
!
        CALL parcom2_seg(coord_g(1,3),
     &                   coord_g(1,4),
     &                   coord_g(1,4), ! NO EFFECT FOR THIS ONE
     &                   nseg,1,1,2,mesh,1,11)

      ENDIF
!
!     SECOND PART:
!     RETRIEVE CMI: CMI(1,ISEG) AND CMI(2,ISEG) ARE X AND Y OF THE
!                   MID-POINT G1G2
!     RETRIEVE JMI: JMI(ISEG) IS THE NUMBER IF TRIANGLE TO WHICH
!                   BELONGS CMI
!
      DO ielem=1,nelem
        DO i = 1,3
          iseg = eltseg(ielem,i)
          nb1 = gloseg(iseg,1)
          nb2 = gloseg(iseg,2)
          xg1 = coord_g(iseg,1)
          yg1 = coord_g(iseg,2)
          xg2 = coord_g(iseg,3)
          yg2 = coord_g(iseg,4)
          ! CENTER OF SEGMENT G1G2
          x_midpoint = 0.5d0*(xg1+xg2)
          y_midpoint = 0.5d0*(yg1+yg2)
          cmi(1,iseg) = x_midpoint
          cmi(2,iseg) = y_midpoint
          ! JMI: IN WHICH ELEMENT IS CMI
          ! WE USE CROSS PRODUCT
          xsom(1)=x(nb1)
          xsom(2)=x(nb2)
          ysom(1)=y(nb1)
          ysom(2)=y(nb2)

          det = (xsom(2)-xsom(1))*(y_midpoint-ysom(1))
     &         -(ysom(2)-ysom(1))*(x_midpoint-xsom(1))

          IF(oriseg(ielem,i).EQ.2) det = -det
!         BE CARREFUL FIXME:
!           - IN CASE OF ISEG IS AT THE INTERFACE OF THE TRIANGLE
!             EACH SUBDOMAIN WILL SEE IT
          IF(det.GE.0.d0) THEN
            jmi(iseg) = ielem
          ENDIF
        ENDDO
      ENDDO
!
!     FOR PARALLELISM
!
      IF(ncsize.GT.1) THEN
!       TODO: Manually creating temporary array
        ALLOCATE(tmp_cmi1(nseg))
        ALLOCATE(tmp_cmi2(nseg))
        tmp_cmi1 = cmi(1,1:nseg)
        tmp_cmi2 = cmi(2,1:nseg)
        CALL parcom2_seg(tmp_cmi1,
     &                   tmp_cmi2,
     &                   cmi,
     &                   nseg,1,1,2,mesh,1,11)
        cmi(1,1:nseg) = tmp_cmi1
        cmi(2,1:nseg) = tmp_cmi2
        DEALLOCATE(tmp_cmi1)
        DEALLOCATE(tmp_cmi2)
      ENDIF
!
! BOUNDARY SEGMENTS
!
      DO ielem=1,nelem
        i1 = ikle(ielem,1)
        i2 = ikle(ielem,2)
        i3 = ikle(ielem,3)
        DO i = 1,3
          iseg = eltseg(ielem,i)
          IF(ifabor(ielem,i).EQ.-1.OR. ifabor(ielem,i).EQ. 0) THEN
!           BOUNDARY SEGMENT
            nb1=gloseg(iseg,1)
            nb2=gloseg(iseg,2)
!           CMI
            cmi(1,iseg)=0.5d0*(x(nb1)+x(nb2))
            cmi(2,iseg)=0.5d0*(y(nb1)+y(nb2))
!           JMI
            jmi(iseg)=ielem
!           ISEG HAS BEEN TREATED
          ENDIF
        ENDDO
      ENDDO
!
! TO VERIFY THAT IT IS WELL DONE
!
      IF(ncsize.LE.1)THEN
        DO iseg=1,nseg
          IF(jmi(iseg).EQ.0)THEN
            WRITE(lu,*)'CENTRE_MASS_SEG: PROBLEM JMI NOT GOOD'
            WRITE(lu,*)'FOR SEGMENT :',iseg
            WRITE(lu,*)'PTS SEG',gloseg(iseg,1),gloseg(iseg,2)
            CALL plante(1)
            stop
          ENDIF
        ENDDO
      ENDIF
!
!---------------------------------------------------------------------
!
      RETURN
      END