rpi_invr.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\tomawac\rpi_invr.f
00002 !
00046                         SUBROUTINE RPI_INVR
00047 !                       *******************
00048 !
00049      &( X     , Y     , NEIGB , NB_CLOSE, RK_D , RX_D , RY_D  , RXX_D ,
00050      &  RYY_D , NPOIN2, I     , QUO   , AC    , MAXNSP, MINDIST )
00051 !
00052 !***********************************************************************
00053 ! TOMAWAC   V6P2                                   25/06/2012
00054 !***********************************************************************
00055 !
00056 !
00057 !
00058 !
00059 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00060 !| AC             |-->| CONSTANT FOR RADIAL FUNCTION COMPUT.
00061 !| I              |-->| POINT INDEX
00062 !| MAXNSP         |-->| CONSTANT FOR MESHFREE TECHNIQUE
00063 !| MINDIST        |-->| CONSTANT FOR RADIAL FUNCTION COMPUT.
00064 !| NB_CLOSE       |-->| ARRAY USED IN THE MESHFREE TECHNIQUE
00065 !| NEIGB          |-->| NEIGHBOUR POINTS FOR MESHFREE METHOD
00066 !| NPOIN2         |-->| NUMBER OF POINTS IN 2D MESH
00067 !| QUO            |-->| CONSTANT FOR RADIAL FUNCTION COMPUT.
00068 !| RK_D           |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
00069 !| RX_D           |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
00070 !| RXX_D          |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
00071 !| RY_D           |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
00072 !| RYY_D          |<->| ARRAY USED IN THE MESHFREE TECHNIQUE
00073 !| X              |-->| ABSCISSAE OF POINTS IN THE MESH
00074 !| Y              |-->| ORDINATES OF POINTS IN THE MESH
00075 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00076 !
00077       IMPLICIT NONE
00078 !
00079 !     VARIABLES IN ARGUMENT
00080 !
00081       INTEGER NPOIN2, MAXNSP, I
00082       INTEGER NEIGB(NPOIN2,MAXNSP),NB_CLOSE(NPOIN2)
00083 !
00084       DOUBLE PRECISION QUO, AC
00085       DOUBLE PRECISION X(NPOIN2), Y(NPOIN2)
00086       DOUBLE PRECISION MINDIST(NPOIN2)
00087       DOUBLE PRECISION RK_D(MAXNSP)
00088       DOUBLE PRECISION RX_D(MAXNSP), RY_D(MAXNSP)
00089       DOUBLE PRECISION RXX_D(MAXNSP), RYY_D(MAXNSP)
00090 !
00091 !     LOCAL VARIABLES
00092 !
00093       INTEGER IP, IPOIN, IP1, IPOIN1, NP
00094 
00095       DOUBLE PRECISION, ALLOCATABLE:: RK_I(:,:), RN(:,:)
00096       DOUBLE PRECISION, ALLOCATABLE:: RX_I(:,:), RY_I(:,:)
00097       DOUBLE PRECISION, ALLOCATABLE:: RXX_I(:,:), RYY_I(:,:)
00098       DOUBLE PRECISION, ALLOCATABLE:: RAD1(:,:)
00099 !
00100       DOUBLE PRECISION DC, WZ, WZX1, WZY1, WZX2, WZY2
00101 !
00102       LOGICAL DEJA
00103       DATA DEJA /.FALSE./
00104 !
00105       SAVE
00106 !
00107 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00108 !
00109       IF(.NOT.DEJA)THEN
00110         ALLOCATE(RK_I(MAXNSP,MAXNSP))
00111         ALLOCATE(RN(MAXNSP,MAXNSP))
00112         ALLOCATE(RX_I(MAXNSP,MAXNSP))
00113         ALLOCATE(RY_I(MAXNSP,MAXNSP))
00114         ALLOCATE(RXX_I(MAXNSP,MAXNSP))
00115         ALLOCATE(RYY_I(MAXNSP,MAXNSP))
00116         ALLOCATE(RAD1(MAXNSP,MAXNSP))
00117         DEJA=.TRUE.
00118       ENDIF
00119 !
00120       DO IP1 =1,NB_CLOSE(I)
00121         IP=NEIGB(I,IP1)
00122         DO IPOIN1 =1,NB_CLOSE(I)
00123           IPOIN=NEIGB(I,IPOIN1)
00124           RAD1(IP1,IPOIN1)=(X(IP)-X(IPOIN))**2+(Y(IP)-Y(IPOIN))**2
00125         ENDDO
00126       ENDDO
00127 !
00128       DO IP1 =1,NB_CLOSE(I)
00129         IP=NEIGB(I,IP1)
00130         DC=MINDIST(I)
00131         DO IPOIN1 =1,NB_CLOSE(I)
00132           IPOIN=NEIGB(I,IPOIN1)
00133           RK_I(IP1,IPOIN1)=(RAD1(IP1,IPOIN1)+(AC*DC)**2)**QUO
00134 !
00135 !         First derivative
00136 !
00137           RY_I(IP1,IPOIN1)=2.D0*QUO*(RAD1(IP1,IPOIN1)+
00138      &     (AC*DC)**2)**(QUO-1.D0)*(Y(IP)-Y(IPOIN))
00139 
00140           RX_I(IP1,IPOIN1)=2.*QUO*(RAD1(IP1,IPOIN1)+
00141      &     (AC*DC)**2)**(QUO-1.D0)*(X(IP)-X(IPOIN))
00142 !
00143 !         Second derivative
00144 !
00145           RYY_I(IP1,IPOIN1) =2.D0*QUO*(RAD1(IP1,IPOIN1)+
00146      &(AC*DC)**2)**(QUO-1.D0)+4.D0*QUO*(QUO-1.D0)*
00147      &(RAD1(IP1,IPOIN1)+(AC*DC)**2)**(QUO-2.D0)*(Y(IP)-Y(IPOIN))**2
00148 !
00149           RXX_I(IP1,IPOIN1)=2.D0*QUO*(RAD1(IP1,IPOIN1)+
00150      &(AC*DC)**2)**(QUO-1.D0)+4.D0*QUO*(QUO-1.D0)*
00151      &(RAD1(IP1,IPOIN1)+(AC*DC)**2)**(QUO-2.D0)*(X(IP)-X(IPOIN))**2
00152         ENDDO
00153       ENDDO
00154 !
00155       RN=RK_I
00156       NP=NB_CLOSE(I)
00157 !
00158       CALL INVERT(RN,NP,MAXNSP)
00159 !
00160       DO IP1 =1,NB_CLOSE(I)
00161         WZ=0.D0
00162         WZX1=0.D0
00163         WZY1=0.D0
00164         WZX2=0.D0
00165         WZY2=0.D0
00166         DO IPOIN1 =1,NB_CLOSE(I)
00167           WZ=WZ+RK_I(1,IPOIN1)*RN(IPOIN1,IP1)
00168           WZX1=WZX1+RX_I(1,IPOIN1)*RN(IPOIN1,IP1)
00169           WZY1=WZY1+RY_I(1,IPOIN1)*RN(IPOIN1,IP1)
00170           WZX2=WZX2+RXX_I(1,IPOIN1)*RN(IPOIN1,IP1)
00171           WZY2=WZY2+RYY_I(1,IPOIN1)*RN(IPOIN1,IP1)
00172         ENDDO
00173 !       write RK etc for the right form for each domain in one row
00174         RK_D(IP1)  = WZ
00175         RX_D(IP1)  = WZX1
00176         RY_D(IP1)  = WZY1
00177         RXX_D(IP1) = WZX2
00178         RYY_D(IP1) = WZY2
00179       ENDDO
00180 !
00181       RETURN
00182       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0