proxim.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE proxim
!                   *****************
!
     &(ip,xp,yp,x,y,np,npoin,ikle,nelem,nelmax)
!
!***********************************************************************
! BIEF   V7P3
!***********************************************************************
!
!brief    IDENTIFIES THE POINTS OF THE MESH CLOSEST TO A SET
!+                OF GIVEN POINTS.
!
!history  J-M HERVOUET (LNHE)
!+        03/07/2009
!+        V6P0
!+
!
!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 (LNHE)
!+        16/11/2012
!+        V6P3
!+   Write statements added.
!
!history  J-M HERVOUET (LNHE)
!+        24/12/2013
!+        V7P0
!+   In parallel the distance between given points and points in the
!+   mesh must be now less than 1.D-8 instead of 1.D-4 (to avoid
!+   problems with subroutine ecrspe in Tomawac)
!
!history  J-M HERVOUET (LNHE)
!+        14/11/2014
!+        V7P0
!+   Checking that a point belongs to at least one sub-domain in
!+   parallel mode.
!
!history  J-M HERVOUET (LNHE)
!+        16/06/2016
!+        V7P2
!+   Reverting to the original algorithm in sequential mode, for parallel
!+   also: choosing the nearest point in the mesh. In case of parallelism
!+   the nearest of all points in all sub-domains is chosen.
!
!history  J-M HERVOUET (LNHE)
!+        29/10/2017
!+        V7P3
!+   Once a processor has found a source point, it must be communicated
!+   to others in case this point is on an interface. This is done by
!+   simply comparing the coordinates.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IKLE           |-->| CONNECTIVITY TABLE.
!| IP             |<--| ADDRESSES OF NEAREST POINTS
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NP             |-->| NUMBER OF POINTS IN THE SET
!| NPOIN          |-->| NUMBER OF POINTS IN THE MESH
!| X              |-->| ABSCISSAE OF POINTS IN THE MESH
!| XP             |-->| ABSCISSAE OF POINTS IN THE SET
!| Y              |-->| ORDINATES OF POINTS IN THE MESH
!| YP             |-->| ORDINATES OF POINTS IN THE SET
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_parallel
      USE bief, ex_proxim => proxim
!
      USE declarations_special
      IMPLICIT NONE
!
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)    :: NP,NPOIN,NELEM,NELMAX
      INTEGER, INTENT(INOUT) :: IP(np)
      INTEGER, INTENT(IN)    :: IKLE(nelmax,3)
!
      DOUBLE PRECISION, INTENT(IN) :: XP(np),YP(np),X(npoin),Y(npoin)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,K,IELEM
      DOUBLE PRECISION X1,Y1,X2,Y2,X3,Y3,A31,A12,A23,DIST2,D2
      DOUBLE PRECISION XX,YY
!
      INTRINSIC sqrt
!
!-----------------------------------------------------------------------
!
      IF(np.GT.0) THEN
      DO k=1,np
        ip(k)=0
        dist2=1.d10
        xx=-1.d10
        yy=-1.d10
!
!       LOOP ON THE TRIANGLES TO FIND THE NEAREST POINT IN THE (SUB-)DOMAIN:
!
        DO ielem=1,nelem
          x1=x(ikle(ielem,1))
          x2=x(ikle(ielem,2))
          x3=x(ikle(ielem,3))
          y1=y(ikle(ielem,1))
          y2=y(ikle(ielem,2))
          y3=y(ikle(ielem,3))
          a31=xp(k)*y3-yp(k)*x3+x3*y1-x1*y3+x1*yp(k)-xp(k)*y1
          a12=xp(k)*y1-yp(k)*x1+x1*y2-x2*y1+x2*yp(k)-xp(k)*y2
          a23=xp(k)*y2-yp(k)*x2+x2*y3-x3*y2+x3*yp(k)-xp(k)*y3
          IF(a31.GT.-1.d-6.AND.a12.GT.-1.d-6.AND.a23.GT.-1.d-6) THEN
!           TAKES THE NEAREST NODE
            DO i=1,3
              d2=(xp(k)-x(ikle(ielem,i)))**2+(yp(k)-y(ikle(ielem,i)))**2
              IF(d2.LT.dist2) THEN
                ip(k)=ikle(ielem,i)
                xx=x(ikle(ielem,i))
                yy=y(ikle(ielem,i))
                dist2=d2
              ENDIF
            ENDDO
          ENDIF
        ENDDO
!
!       CHECKING THAT THE POINT IS IN THE GLOBAL DOMAIN
!       IF YES PRINTING THE COORDINATES OF THE NEAREST POINT FOUND
!
        i=ip(k)
        IF(ncsize.GT.1) i=p_max(i)
        IF(i.EQ.0) THEN
!         THE POINT IS NOT IN THE DOMAIN
          WRITE(lu,*)
          WRITE(lu,*) 'SPECTRUM OR SOURCE POINT ',k,' OUTSIDE DOMAIN'
          CALL plante(1)
          stop
        ELSE
!         THE POINT IS IN THE DOMAIN. IN PARALLEL SEVERAL NEAREST POINTS
!         MAY HAVE BEEN FOUND, FINDING THE REAL NEAREST ONE
          x1=xx
          y1=yy
          IF(ncsize.GT.1) THEN
            IF(dist2.NE.p_min(dist2)) THEN
              x1=0.d0
              y1=0.d0
            ENDIF
            x1=p_min(x1)+p_max(x1)
            y1=p_min(y1)+p_max(y1)
!           ALL PROCESSORS HAVING THE POINT MUST KNOW IT
!           YET POINTS WITHIN A TRIANGLE IN ANOTHER PROCESSOR
!           AND FALLEN BACK ON AN INTERFACE HAVE BEEN OVERLOOKED
!           SO FAR BY OTHER PROCESSORS.
            DO i=1,npoin
              IF(x(i).EQ.x1.AND.y(i).EQ.y1) ip(k)=i
            ENDDO
          ENDIF
          WRITE(lu,*)
          WRITE(lu,*) 'SOURCE POINT ',k,'PUT ON POINT'
          WRITE(lu,*) x1,' AND ',y1
          d2 = sqrt(p_min(dist2))
          WRITE(lu,*) 'LOCATED AT ',d2,' METRES'
!         LINE FEED FOR THE LISTING
          IF(k.EQ.np) WRITE(lu,*)
        ENDIF
      ENDDO
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END