checkmesh.f

Go to the documentation of this file.

!                   ********************
                    SUBROUTINE checkmesh
!                   ********************
!
     &(mesh,npoin)
!
!***********************************************************************
! BIEF   V7P1
!***********************************************************************
!
!brief    Checks the mesh.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        23/10/2015
!+        V7P1
!+   First version. Checking the coordinates and the number of
!+   neighbours of points.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| MESH           |-->| MESH STRUCTURE.
!| NPOIN          |-->| NUMBER OF POINTS (NPOIN2 IN 3D)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_checkmesh => checkmesh
      USE interface_parallel
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          , INTENT(IN)    :: NPOIN
      TYPE(bief_mesh)  , INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,J,IDISTMIN,JDISTMIN,IELEM,I1,I2,I3,OVER
      DOUBLE PRECISION, POINTER :: X(:),Y(:)
      DOUBLE PRECISION DIST2,DIST2MIN
      INTRINSIC sqrt,nint
      LOGICAL STOP_ERROR
!
!-----------------------------------------------------------------------
!
      x=>mesh%X%R
      y=>mesh%Y%R
      stop_error=.false.
!
!-----------------------------------------------------------------------
!     CHECKING THE COORDINATES
!-----------------------------------------------------------------------
!
      dist2min=1.d20
      DO i=1,npoin-1
        DO j=i+1,npoin
          dist2=(x(i)-x(j))**2+(y(i)-y(j))**2
          IF(dist2.LT.dist2min) THEN
            dist2min=dist2
            idistmin=i
            jdistmin=j
          ENDIF
        ENDDO
      ENDDO
      IF(ncsize.GT.1) THEN
        IF(p_min(dist2min).EQ.dist2min) THEN
          idistmin=mesh%KNOLG%I(idistmin)
          jdistmin=mesh%KNOLG%I(jdistmin)
        ELSE
          idistmin=0
          jdistmin=0
          dist2min=0.d0
        ENDIF
        dist2min=p_max(dist2min)
        idistmin=p_max(idistmin)
        jdistmin=p_max(jdistmin)
      ENDIF
      WRITE(lu,*)
      WRITE(lu,*) 'CHECKING THE MESH'
      WRITE(lu,*)
      WRITE(lu,*) 'SMALLEST DISTANCE BETWEEN TWO POINTS:',
     &            sqrt(dist2min)
      WRITE(lu,*) 'BETWEEN POINTS: ',idistmin,' AND ',jdistmin
      IF(dist2min.LT.1.d-3) THEN
        WRITE(lu,*) 'VALUE TOO SMALL'
        stop_error=.true.
      ENDIF
!
!-----------------------------------------------------------------------
!     CHECKING THE NUMBER OF NEIGHBOURS
!-----------------------------------------------------------------------
!
      DO i=1,npoin
        mesh%T%R(i)=0.d0
      ENDDO
!
      DO ielem=1,mesh%NELEM
        i1=mesh%IKLE%I(ielem)
        mesh%T%R(i1)=mesh%T%R(i1)+2.d0
        i2=mesh%IKLE%I(ielem+mesh%NELMAX)
        mesh%T%R(i2)=mesh%T%R(i2)+2.d0
        i3=mesh%IKLE%I(ielem+2*mesh%NELMAX)
        mesh%T%R(i3)=mesh%T%R(i3)+2.d0
      ENDDO
      IF(ncsize.GT.1) CALL parcom(mesh%T,2,mesh)
!
      over=0
!
      DO i=1,npoin
        IF(mesh%T%R(i).LT.3.d0) THEN
          IF(ncsize.GT.1) THEN
            j=mesh%KNOLG%I(i)
          ELSE
            j=i
          ENDIF
          WRITE(lu,*) 'POINT ',j,' HAS ONLY ',
     &                nint(mesh%T%R(i)),' NEIGHBOURS'
          IF(mesh%T%R(i).LT.1.d0) THEN
            stop_error=.true.
          ELSE
            over=over+1
          ENDIF
        ENDIF
      ENDDO
!
!     THIS SUM IS CORRECT BECAUSE A POINt WITH 2 NEIGHBOURS CANNOT
!     BELONG TO 2 SUB-DOMAINS
      IF(ncsize.GT.1) over=p_sum(over)
!
      IF(over.GT.1) THEN
        WRITE(lu,*) over,' OVERCONSTRAINED TRIANGLES'
        WRITE(lu,*) 'POSSIBLE BUT NOT RECOMMENDED'
      ENDIF
!
!-----------------------------------------------------------------------
!     STOPPING IF FATAL ERROR
!-----------------------------------------------------------------------
!
      IF(stop_error) THEN
        WRITE(lu,*) 'ISSUE DETECTED WHILE CHECING THE MESH'
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END