centre_mass_seg.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\centre_mass_seg.f
00002 !
00061                         SUBROUTINE CENTRE_MASS_SEG
00062 !                       **************************
00063 !
00064      &(X,Y,COORD_G,IKLE,NPOIN,ELTSEG,ORISEG,NELEM,NSEG,JMI,CMI,GLOSEG,
00065      & IFABOR,MESH)
00066 !
00067 !***********************************************************************
00068 ! BIEF   V6P3                                                18/12/2012
00069 !***********************************************************************
00070 !
00071 !           LEFT OF SEGMENTS I.E. FOR A SEGMENT ISEG:
00072 !            - COORD_G(1,ISEG) AND COORD_G(2,ISEG) ARE THE COORDINATES OF
00073 !              THE CENTRE OF GRAVITY OF THE FIRST NEIGHBORING TRIANGLE
00074 !            - COORD_G(3,ISEG) AND COORD_G(4,ISEG) ARE THE COORDINATES OF
00075 !              THE CENTRE OF GRAVITY OF THE SECOND NEIGHBORING TRIANGLE
00076 !         GIVES COORDINATES OF THE MIDDLE POINT OF SEGMENT G1G2 (STOCKED
00077 !             AT CMI(1,ISEG) AND CMI(2,ISEG)) AND THE ELEMENT NUMBER TO
00078 !             TO WHICH BELONGS CMI
00079 !
00080 !
00081 !
00082 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00083 !| CMI            |<--| COORDINATES OF MID-INTERFACE POINTS
00084 !| COORD_G        |<--| CENTER OF MASS OF ELEMENTS NEIGHBORS OF AN EDGE
00085 !|                |   |  COORD_G(1,ISEG) AND COORD_G(2,ISEG) ARE X1 AND Y1
00086 !|                |   |  COORD_G(3,ISEG) AND COORD_G(4,ISEG) ARE X2 AND Y2
00087 !| ELTSEG         |-->| SEGMENTS FORMING AN ELEMENT
00088 !| GLOSEG         |-->|  GLOBAL NUMBERS OF VERTICES OF SEGMENTS
00089 !| IKLE           |-->| CONNECTIVITY TABLE.
00090 !| IFABOR         |-->| IFABOR(IEL,I) IS THE ELEMENT BEHIND THE EDGE I OF
00091 !|                |---|  ELEMENT IEL, OTHERWISE
00092 !|                |---|     IFABOR(IEL,I) = -2 : THIS IS INTERFACE EDGE
00093 !|                |---|     IFABOR(IEL,I) = 0  : THIS IS BOUNDARY EDGE
00094 !|                |---|     IFABOR(IEL,I) = -1 : THIS IS LIQUID BOUNDARY EDGE
00095 !| JMI            |<--| NUMBER OF TRIANGLE TO WHICH BELONGS THE
00096 !|                |   |       MID-INTERFACE POINT.
00097 !| MESH           |-->| MESH STRUCTURE
00098 !| NELEM          |-->| NUMBER OF ELEMENTS
00099 !| NPOIN          |-->| NUMBER OF POINTS
00100 !| NSEG           |-->| NUMBER OF SEGMENTS
00101 !|                |---|
00102 !| ORISEG         |-->| ORIENTATION OF SEGMENTS FORMING AN
00103 !|                |   |        ELEMENT 1:ANTI 2:CLOCKWISE
00104 !| X              |-->| ABSCISSAE OF POINTS IN THE MESH
00105 !| Y              |-->| ORDINATES OF POINTS IN THE MESH
00106 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00107 !
00108       USE BIEF, ONLY : INPOLY
00109       USE BIEF_DEF
00110       IMPLICIT  NONE
00111       INTEGER LNG,LU
00112       COMMON/INFO/LNG,LU
00113 !
00114 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00115 !
00116       INTEGER, INTENT(IN)             :: NSEG,NELEM,NPOIN
00117       INTEGER, INTENT(IN)             :: IKLE(NELEM,3)
00118       INTEGER, INTENT(IN)             :: ELTSEG(NELEM,3)
00119       INTEGER, INTENT(IN)             :: ORISEG(NELEM,3)
00120       INTEGER, INTENT(INOUT)          :: JMI(NSEG)
00121       INTEGER, INTENT(IN)             :: GLOSEG(NSEG,2)
00122       DOUBLE PRECISION, INTENT(INOUT) :: CMI(2,NSEG)
00123       DOUBLE PRECISION, INTENT(INOUT) :: COORD_G(NSEG,4)
00124       DOUBLE PRECISION, INTENT(IN)    :: X(NPOIN),Y(NPOIN)
00125       INTEGER, INTENT(IN)             :: IFABOR(NELEM,3)
00126       TYPE(BIEF_MESH), INTENT(INOUT)  :: MESH
00127 !
00128 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00129 !
00130       INTEGER ISEG,NB1,NB2,IELEM,I,I1,I2,I3
00131       DOUBLE PRECISION XEL,YEL,XG1,XG2,YG1,YG2,XSOM(3),YSOM(3)
00132       DOUBLE PRECISION X_MIDPOINT,Y_MIDPOINT
00133       LOGICAL DEJA
00134 !
00135 !-----------------------------------------------------------------------
00136 !   INITIALIZATION OF VARIABLES
00137 !
00138       DO I=1,NSEG
00139         COORD_G(I,1) = 0.D0
00140         COORD_G(I,2) = 0.D0
00141         COORD_G(I,3) = 0.D0
00142         COORD_G(I,4) = 0.D0
00143         CMI(1,I)     = 0.D0
00144         CMI(2,I)     = 0.D0
00145         JMI(I)       = 0
00146       ENDDO
00147       IF(NCSIZE.GT.1)THEN
00148         DO I=1,NSEG
00149           MESH%MSEG%X%R(I)=0.D0
00150         ENDDO
00151       ENDIF
00152 !    RECUPERATE THE COORDINATES OF CENTER OF GRAVITY OF THE TRIANGLES
00153 !    AT THE RIGHT AND AT THE LEFT OF THE SEGMENT BETWEEN TWO NODES
00154 !
00155       DO IELEM=1, NELEM
00156         I1 = IKLE(IELEM,1)
00157         I2 = IKLE(IELEM,2)
00158         I3 = IKLE(IELEM,3)
00159         XEL = (X(I1)+X(I2)+X(I3))/3.0D0
00160         YEL = (Y(I1)+Y(I2)+Y(I3))/3.0D0
00161         DO I = 1,3
00162           ISEG = ELTSEG(IELEM,I)
00163           IF(ORISEG(IELEM,I).EQ.1)THEN ! G IS LEFT OF THE EDGE
00164             COORD_G(ISEG,1) =XEL
00165             COORD_G(ISEG,2) =YEL
00166           ELSEIF(ORISEG(IELEM,I).EQ.2)THEN !G IS RIGHT OF THE EDGE
00167             COORD_G(ISEG,3) =XEL
00168             COORD_G(ISEG,4) =YEL
00169           ENDIF
00170         ENDDO
00171       ENDDO
00172 !
00173       IF(NCSIZE.GT.1) THEN
00174         CALL PARCOM2_SEG(COORD_G(1,1),
00175      &                   COORD_G(1,2),
00176      &                   COORD_G(1,2), ! NO EFFECT FOR THIS ONE
00177      &                   NSEG,1,1,2,MESH,1,11)
00178 !
00179         CALL PARCOM2_SEG(COORD_G(1,3),
00180      &                   COORD_G(1,4),
00181      &                   COORD_G(1,4), ! NO EFFECT FOR THIS ONE
00182      &                   NSEG,1,1,2,MESH,1,11)
00183 
00184       ENDIF
00185 !
00186 !     SECOND PART:
00187 !     RETRIEVE CMI: CMI(1,ISEG) AND CMI(2,ISEG) ARE X AND Y OF THE
00188 !                   MID-POINT G1G2
00189 !     RETRIEVE JMI: JMI(ISEG) IS THE NUMBER IF TRIANGLE TO WHICH
00190 !                   BELONGS CMI
00191 !
00192       DO IELEM=1,NELEM
00193         I1 = IKLE(IELEM,1)
00194         I2 = IKLE(IELEM,2)
00195         I3 = IKLE(IELEM,3)
00196         DO I = 1,3
00197           ISEG = ELTSEG(IELEM,I)
00198           NB1 = GLOSEG(ISEG,1)
00199           NB2 = GLOSEG(ISEG,2)
00200           XG1 = COORD_G(ISEG,1)
00201           YG1 = COORD_G(ISEG,2)
00202           XG2 = COORD_G(ISEG,3)
00203           YG2 = COORD_G(ISEG,4)
00204           ! CENTER OF SEGMENT G1G2
00205           X_MIDPOINT = 0.5D0*(XG1+XG2)
00206           Y_MIDPOINT = 0.5D0*(YG1+YG2)
00207           ! JMI: IN WHICH ELEMENT IS CMI
00208           ! WE USE INPOLY FUNCTION (SEE INPOLY.F)
00209           XSOM(1)=X(I1)
00210           XSOM(2)=X(I2)
00211           XSOM(3)=X(I3)
00212           YSOM(1)=Y(I1)
00213           YSOM(2)=Y(I2)
00214           YSOM(3)=Y(I3)
00215           IF(INPOLY(X_MIDPOINT,Y_MIDPOINT,XSOM,YSOM,3 ))THEN
00216             CMI(1,ISEG) = X_MIDPOINT
00217             CMI(2,ISEG) = Y_MIDPOINT
00218             JMI(ISEG)= IELEM
00219           ENDIF
00220         ENDDO
00221       ENDDO
00222 !
00223 !     FOR PARALLELISM
00224 !
00225       IF(NCSIZE.GT.1) THEN
00226 !       NOTE JMH: CMI(1,1:NSEG) HERE IMPLIES A TEMPORARY COPY BY THE COMPILER
00227 !                 AND THEN A REDISTRIBUTION
00228         CALL PARCOM2_SEG(CMI(1,1:NSEG),
00229      &                   CMI(2,1:NSEG),
00230      &                   CMI,
00231      &                   NSEG,1,1,2,MESH,1,11)
00232       ENDIF
00233 !
00234 !     TEST TO SEE IF ALL INTERNAL EDGES ARE WELL TREATED
00235 !     NOTE RA: FOR MALPASSET SMALL, MESH IS VERY POOR AND WITHOUT
00236 !     THE FOLLOWING TEST, IT CRASHES
00237       DEJA =.FALSE.
00238       DO IELEM=1,NELEM
00239         I1 = IKLE(IELEM,1)
00240         I2 = IKLE(IELEM,2)
00241         I3 = IKLE(IELEM,3)
00242         DO I = 1,3
00243           ISEG = ELTSEG(IELEM,I)
00244           IF(JMI(ISEG).EQ.0.AND.
00245      &       IFABOR(IELEM,I).NE.-1.AND.IFABOR(IELEM,I).NE.0) THEN
00246             IF(.NOT.DEJA)THEN
00247               WRITE(LU,*)'MESH WITH POOR QUALITY '
00248               WRITE(LU,*)'FOR INSTANCE, SEE ELEMENT :',IELEM
00249               WRITE(LU,*)'WITH NODES',I1,I2,I3
00250               DEJA=.TRUE.
00251             ENDIF
00252           ENDIF
00253           CMI(1,ISEG)=0.5D0*(COORD_G(ISEG,1)+COORD_G(ISEG,3))
00254           CMI(2,ISEG)=0.5D0*(COORD_G(ISEG,2)+COORD_G(ISEG,4))
00255           JMI(ISEG)=IELEM
00256         ENDDO
00257       ENDDO
00258 !
00259 ! BOUNDARY SEGMENTS
00260 !
00261       DO IELEM=1,NELEM
00262         I1 = IKLE(IELEM,1)
00263         I2 = IKLE(IELEM,2)
00264         I3 = IKLE(IELEM,3)
00265         DO I = 1,3
00266           ISEG = ELTSEG(IELEM,I)
00267           IF(IFABOR(IELEM,I).EQ.-1.OR. IFABOR(IELEM,I).EQ. 0) THEN
00268 !           BOUNDARY SEGMENT
00269             NB1=GLOSEG(ISEG,1)
00270             NB2=GLOSEG(ISEG,2)
00271 !           CMI
00272             CMI(1,ISEG)=0.5D0*(X(NB1)+X(NB2))
00273             CMI(2,ISEG)=0.5D0*(Y(NB1)+Y(NB2))
00274 !           JMI
00275             JMI(ISEG)=IELEM
00276 !           ISEG HAS BEEN TREATED
00277           ENDIF
00278         ENDDO
00279       ENDDO
00280 !
00281 ! TO VERIFY THAT IT IS WELL DONE
00282 !
00283       IF(NCSIZE.LE.1)THEN
00284         DO ISEG=1,NSEG
00285           IF(JMI(ISEG).EQ.0)THEN
00286             WRITE(LU,*)'CENTRE_MASS_SEG: PROBLEM JMI NOT GOOD'
00287             WRITE(LU,*)'FOR SEGMENT :',ISEG
00288             CALL PLANTE(1)
00289             STOP
00290           ENDIF
00291         ENDDO
00292       ENDIF
00293 !
00294 !---------------------------------------------------------------------
00295 !
00296       RETURN
00297       END

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