rpi_invr.f

Go to the documentation of this file.

!                       *******************
                        SUBROUTINE rpi_invr
!                       *******************
!
     &( x     , y     , neigb , nb_close, rk_d , rx_d , ry_d  , rxx_d ,
     &  ryy_d , npoin2, i     , quo   , ac    , maxnsp, mindist )
!
!***********************************************************************
! TOMAWAC   V6P2                                   25/06/2012
!***********************************************************************
!
!brief    DIFFRACTION
!+         CALCULATION OF THE RADIAL FUNCTION FOR THE
!+         FREE-MESH METHOD
!
!history  E. KRIEZI (LNH)
!+        04/12/2006
!+        V5P5
!+
!
!history  G.MATTAROLO (EDF - LNHE)
!+        23/06/2012
!+        V6P2
!+   Modification for V6P2
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AC             |-->| CONSTANT FOR RADIAL FUNCTION COMPUT.
!| I              |-->| POINT INDEX
!| MAXNSP         |-->| CONSTANT FOR MESHFREE TECHNIQUE
!| MINDIST        |-->| CONSTANT FOR RADIAL FUNCTION COMPUT.
!| NB_CLOSE       |-->| ARRAY USED IN THE MESHFREE TECHNIQUE
!| NEIGB          |-->| NEIGHBOUR POINTS FOR MESHFREE METHOD
!| NPOIN2         |-->| NUMBER OF POINTS IN 2D MESH
!| QUO            |-->| CONSTANT FOR RADIAL FUNCTION COMPUT.
!| RK_D           |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
!| RX_D           |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
!| RXX_D          |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
!| RY_D           |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
!| RYY_D          |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
!| X              |-->| ABSCISSAE OF POINTS IN THE MESH
!| Y              |-->| ORDINATES OF POINTS IN THE MESH
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_tomawac, ONLY : deja_rpi,rk_i,rn,rx_i,ry_i,
     &                                 rxx_i,ryy_i,rad1
      USE interface_tomawac, ex_rpi_invr => rpi_invr
      IMPLICIT NONE
!
!     VARIABLES IN ARGUMENT
!
      INTEGER, INTENT(IN)    :: NPOIN2, MAXNSP, I
      INTEGER, INTENT(IN)    :: NEIGB(npoin2,maxnsp), NB_CLOSE(npoin2)
      DOUBLE PRECISION, INTENT(IN)    :: QUO, AC
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin2), Y(npoin2)
      DOUBLE PRECISION, INTENT(IN)    :: MINDIST(npoin2)
      DOUBLE PRECISION, INTENT(INOUT) :: RK_D(maxnsp)
      DOUBLE PRECISION, INTENT(INOUT) :: RX_D(maxnsp), RY_D(maxnsp)
      DOUBLE PRECISION, INTENT(INOUT) :: RXX_D(maxnsp), RYY_D(maxnsp)
!
!     LOCAL VARIABLES
!
      INTEGER IP, IPOIN, IP1, IPOIN1, NP

!
      DOUBLE PRECISION DC, WZ, WZX1, WZY1, WZX2, WZY2
!      
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
! 
      IF(.NOT.deja_rpi)THEN
        ALLOCATE(rk_i(maxnsp,maxnsp))
        ALLOCATE(rn(maxnsp,maxnsp))
        ALLOCATE(rx_i(maxnsp,maxnsp))
        ALLOCATE(ry_i(maxnsp,maxnsp))
        ALLOCATE(rxx_i(maxnsp,maxnsp))
        ALLOCATE(ryy_i(maxnsp,maxnsp))
        ALLOCATE(rad1(maxnsp,maxnsp))
        deja_rpi=.true.
      ENDIF
!
      DO ip1 =1,nb_close(i)
        ip=neigb(i,ip1)
        DO ipoin1 =1,nb_close(i)
          ipoin=neigb(i,ipoin1)
          rad1(ip1,ipoin1)=(x(ip)-x(ipoin))**2+(y(ip)-y(ipoin))**2
        ENDDO
      ENDDO
!
      DO ip1 =1,nb_close(i)
        ip=neigb(i,ip1)
        dc=mindist(i)
        DO ipoin1 =1,nb_close(i)
          ipoin=neigb(i,ipoin1)
          rk_i(ip1,ipoin1)=(rad1(ip1,ipoin1)+(ac*dc)**2)**quo
!
!         First derivative
!
          ry_i(ip1,ipoin1)=2.d0*quo*(rad1(ip1,ipoin1)+
     &     (ac*dc)**2)**(quo-1.d0)*(y(ip)-y(ipoin))

          rx_i(ip1,ipoin1)=2.*quo*(rad1(ip1,ipoin1)+
     &     (ac*dc)**2)**(quo-1.d0)*(x(ip)-x(ipoin))
!
!         Second derivative
!
          ryy_i(ip1,ipoin1) =2.d0*quo*(rad1(ip1,ipoin1)+
     &(ac*dc)**2)**(quo-1.d0)+4.d0*quo*(quo-1.d0)*
     &(rad1(ip1,ipoin1)+(ac*dc)**2)**(quo-2.d0)*(y(ip)-y(ipoin))**2
!
          rxx_i(ip1,ipoin1)=2.d0*quo*(rad1(ip1,ipoin1)+
     &(ac*dc)**2)**(quo-1.d0)+4.d0*quo*(quo-1.d0)*
     &(rad1(ip1,ipoin1)+(ac*dc)**2)**(quo-2.d0)*(x(ip)-x(ipoin))**2
        ENDDO
      ENDDO
!
      rn=rk_i
      np=nb_close(i)
!
      CALL invert(rn,np,maxnsp)
!
      DO ip1 =1,nb_close(i)
        wz=0.d0
        wzx1=0.d0
        wzy1=0.d0
        wzx2=0.d0
        wzy2=0.d0
        DO ipoin1 =1,nb_close(i)
          wz=wz+rk_i(1,ipoin1)*rn(ipoin1,ip1)
          wzx1=wzx1+rx_i(1,ipoin1)*rn(ipoin1,ip1)
          wzy1=wzy1+ry_i(1,ipoin1)*rn(ipoin1,ip1)
          wzx2=wzx2+rxx_i(1,ipoin1)*rn(ipoin1,ip1)
          wzy2=wzy2+ryy_i(1,ipoin1)*rn(ipoin1,ip1)
        ENDDO
!       write RK etc for the right form for each domain in one row
        rk_d(ip1)  = wz
        rx_d(ip1)  = wzx1
        ry_d(ip1)  = wzy1
        rxx_d(ip1) = wzx2
        ryy_d(ip1) = wzy2
      ENDDO
!
      RETURN
      END