frmset.f

Go to the documentation of this file.

!                       *****************
                        SUBROUTINE frmset
!                       *****************
!
     &( neigb , npoin2, nelem2,
     &  ikle  , rk    , rx    , ry    , rxx   , ryy )
!
!***********************************************************************
! TOMAWAC   V6P2                                   25/06/2012
!***********************************************************************
!
!brief    DIFFRACTION
!+
!+         SETTING THE DOMAINS FOR THE FREE-MESH METHOD
!
!history  E. KRIEZI (LNH)
!+        04/12/2006
!+        V5P5
!+
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IKLE           |-->| TRANSITION BETWEEN LOCAL AND GLOBAL NUMBERING
!|                |   | OF THE 2D MESH
!| MAXNSP         |-->| CONSTANT FOR MESHFREE TECHNIQUE
!| NEIGB          |<->| NEIGHBOUR POINTS FOR MESHFREE METHOD
!| NELEM2         |-->| NUMBER OF ELEMENTS IN 2D MESH
!| NPOIN2         |-->| NUMBER OF POINTS IN 2D MESH
!| RK             |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
!| RX             |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
!| RXX            |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
!| RY             |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
!| RYY            |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_tomawac, ONLY : ilm_poin,clm,kacc,nb_c,sur_p,
     &                                 stack,stack2,already_pom,
     &                                 mindist, deja_frmset,
     &                                 x, y, nb_close, maxnsp 
      USE declarations_special
      USE interface_tomawac, ex_frmset => frmset
      IMPLICIT NONE
!

!.....VARIABLES IN ARGUMENT
!     """"""""""""""""""""
      INTEGER,INTENT(IN)            :: NPOIN2, NELEM2
      INTEGER,INTENT(INOUT)         :: NEIGB(npoin2,maxnsp)
      INTEGER,INTENT(IN)            :: IKLE(nelem2,3)
      DOUBLE PRECISION,INTENT(INOUT):: RK(maxnsp,npoin2)
      DOUBLE PRECISION,INTENT(INOUT):: RX(maxnsp,npoin2)
      DOUBLE PRECISION,INTENT(INOUT):: RY(maxnsp,npoin2)
      DOUBLE PRECISION,INTENT(INOUT):: RXX(maxnsp,npoin2)
      DOUBLE PRECISION,INTENT(INOUT):: RYY(maxnsp,npoin2)
!
!.....LOCAL VARIABLES
!     """""""""""""""
!| NRD            |-->| CONSTANT FOR MESHFREE TECHNIQUE
      INTEGER, PARAMETER :: NRD = 30

      INTEGER IP, IPOIN, IP2, I
      INTEGER ICLM, J, IELEM,ILM,CLMAX
      INTEGER M, ICST, ICST2, NCST, IP_S, ILP, L(2)
      DOUBLE PRECISION AC,QUO,RAD1
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      IF(.NOT.deja_frmset)THEN
        clmax=0
        DO ipoin =1, npoin2
          iclm=0
          DO  ielem=1,nelem2
            IF (ipoin.EQ.ikle(ielem,1)) THEN
              iclm=iclm+1
            ELSEIF (ipoin.EQ.ikle(ielem,2)) THEN
              iclm=iclm+1
            ELSEIF (ipoin.EQ.ikle(ielem,3)) THEN
              iclm=iclm+1
            ENDIF
          ENDDO
          IF(iclm.GT.clmax) clmax=iclm
        ENDDO
        ALLOCATE(ilm_poin(npoin2,clmax))
        ALLOCATE(clm(npoin2))
        ALLOCATE(kacc(npoin2))
        ALLOCATE(nb_c(npoin2))
        ALLOCATE(sur_p(npoin2,clmax))
        ALLOCATE(stack(npoin2))
        ALLOCATE(stack2(npoin2))
        ALLOCATE(already_pom(npoin2))
        ALLOCATE(mindist(npoin2))
        deja_frmset=.true.
      ENDIF
!
! ILM_POIN array with the elements to which a point belongs
! CLM(IP) array with the number of elements for each point
!  for IP belong to the elements  ILM_POIN(CLM(IP-1)+1:CLM(IP))
!
      DO ipoin =1, npoin2
        iclm=0
        DO  ielem=1,nelem2
          IF (ipoin.EQ.ikle(ielem,1)) THEN
            iclm=iclm+1
            ilm_poin(ipoin,iclm)=ielem
          ELSEIF (ipoin.EQ.ikle(ielem,2)) THEN
            iclm=iclm+1
            ilm_poin(ipoin,iclm)=ielem
          ELSEIF (ipoin.EQ.ikle(ielem,3)) THEN
            iclm=iclm+1
            ilm_poin(ipoin,iclm)=ielem
          ENDIF
        ENDDO
        clm(ipoin)=iclm
        IF(clm(ipoin).GT.clmax) THEN
          WRITE(lu,*)  'POINT',ipoin,'HAS TOO MANY NEIGHBOUR: '
     &          ,clm(ipoin)
          CALL plante(1)
          stop
        ENDIF
      ENDDO
!
! searching for the points which are around the point IPOIN
! and add the to a look up array SUR_P(IPOIN,NB_C(IPOIN))
!
!    Initialize all the arrays  and logics for the new subdomain
      DO ip=1,npoin2
        already_pom(ip) =.false.
      ENDDO
!
      DO ipoin=1,npoin2
        nb_c(ipoin)=0
        mindist(ipoin)=1.e+6
        DO ilm=1,clm(ipoin)
          ielem=ilm_poin(ipoin,ilm)
!         loop over 3 nodes of each triangle
          DO j=1,3
!           test if the selected node belongs to the triangle
            IF (ikle(ielem,j).EQ.ipoin) THEN
              IF (j.EQ.1) THEN
                l(1)=ikle(ielem,2)
                l(2)=ikle(ielem,3)
              ENDIF
              IF (j.EQ.2) THEN
                l(1)=ikle(ielem,1)
                l(2)=ikle(ielem,3)
              ENDIF
              IF (j.EQ.3) THEN
                l(1)=ikle(ielem,1)
                l(2)=ikle(ielem,2)
              ENDIF
            ENDIF
          ENDDO
!
          DO m=1,2
            IF (.NOT.already_pom(l(m))) THEN
              nb_c(ipoin)=nb_c(ipoin)+1
              sur_p(ipoin,nb_c(ipoin)) =l(m)
              already_pom(l(m)) =.true.
            ENDIF
          ENDDO
!
        ENDDO
!
!   CALCULATE DISTANCE of EVERY POINT TO THE NEIGHBOUR POINTS
        DO j=1,nb_c(ipoin)
          ip=sur_p(ipoin,j)
          rad1=sqrt((x(ip)-x(ipoin))**2+(y(ip)-y(ipoin))**2)
          IF(rad1.LE.mindist(ipoin)) mindist(ipoin)=rad1
          already_pom(ip) =.false.
        ENDDO

      ENDDO
!
! make the subdomain search over the nearest point of each point and
!  add them in the NEIGB(IPOIN,MAXNSP) aray
      DO ipoin=1,npoin2
        nb_close(ipoin)=1
        neigb(ipoin,1) =ipoin
        already_pom(ipoin) =.true.
        ncst=1
        stack(ncst)=ipoin
!       WRITE(LU,*) (SUR_P(STACK(NCST),j),j=1,7)
!
! ipoin is the main point of domain ipoin
! around the ipoin do a search in the elements it belongs
! loop around the point of Stack and do what you did before
!
! while loop
        DO
          icst2=0
          DO icst =1,ncst
            ip= stack(icst)

            DO ilp=1,nb_c(ip)
              ip_s=sur_p(ip,ilp)
              IF (.NOT.already_pom(ip_s)) THEN
                nb_close(ipoin)=nb_close(ipoin)+1
                neigb(ipoin,nb_close(ipoin)) =ip_s
                already_pom(ip_s) =.true.
                IF(nb_close(ipoin).GE.nrd) GOTO 222
                icst2=icst2+1
                stack2(icst2)=ip_s
              ENDIF
            ENDDO ! ILP
          ENDDO ! ICST
          ncst=icst2
          stack=stack2
        ENDDO
222   CONTINUE
! end of while loop
!
!subdomain (Ipoin) finish after initializing
! logic goto to the next subdomain
        DO j=1,nb_close(ipoin)! initialize already for points logic
          ip2=neigb(ipoin,j)
          already_pom(ip2) =.false. ! initialize already for points logic
        ENDDO
      ENDDO  !1,NPOIN2
!
! CALCULATE THE RADIAL FUNCTION OF RPI
! AND INVERSE MATRICES OF EACH SUB DOMAIN
!
      quo = 1.03d0
      ac = 8.d0
      DO i=1,npoin2
        CALL rpi_invr(x, y, neigb, nb_close,
     &     rk(1,i), rx(1,i), ry(1,i), rxx(1,i), ryy(1,i),
     &     npoin2, i, quo, ac, maxnsp, mindist)
!
      ENDDO
!
      RETURN
      END