gettriebe.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\gettriebe.f
00002 !
00082                      SUBROUTINE GETTRIEBE
00083 !                    ********************
00084 !
00085      &(XAUX,AD,AX,TETA,IKLE,NPOIN,NELEM,NELMAX,MESH,IELM3,NELEM2,NPLAN,
00086      & KNOLG)
00087 !
00088 !***********************************************************************
00089 ! BIEF   V6P2                                   21/08/2010
00090 !***********************************************************************
00091 !
00092 !
00093 !
00094 !
00095 !
00096 !
00097 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00098 !| AD             |-->| DIAGONAL TERMS OF MATRIX
00099 !| AX             |-->| OFF-DIAGONAL TERMS OF MATRIX
00100 !| IELM3          |-->| TYPE OF ELEMENT
00101 !| IKLE           |-->| CONNECTIVITY TABLE
00102 !| KNOLG          |-->| ORIGINAL (I.E. SCALAR MODE) GLOBAL NUMBER OF POINTS
00103 !| MESH           |-->| MESH STRUCTURE
00104 !| NELEM          |-->| NUMBER OF ELEMENTS
00105 !| NELEM2         |-->| NUMBER OF TRIANGLES OF ORIGINAL 2D MESH
00106 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00107 !| NPLAN          |-->| NUMBER OF PLANES IN THE ORIGINAL MESH
00108 !| NPOIN          |-->| NUMBER OF POINTS
00109 !| TETA           |-->| COEFFICIENT USED IN THE RESULT
00110 !| XAUX           |<--| THE RESULTING MATRIX
00111 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00112 !
00113       USE BIEF, EX_GETTRIEBE => GETTRIEBE
00114       USE DECLARATIONS_TELEMAC, ONLY : TETRA,ISEGT
00115 !
00116       IMPLICIT NONE
00117       INTEGER LNG,LU
00118       COMMON/INFO/LNG,LU
00119 !
00120 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00121 !
00122       INTEGER, INTENT(IN) :: NELEM,NELMAX,NPOIN,IELM3,NELEM2,NPLAN
00123       INTEGER, INTENT(IN) :: IKLE(NELMAX,*),KNOLG(NPOIN)
00124 !
00125       DOUBLE PRECISION, INTENT(IN)    :: TETA
00126       DOUBLE PRECISION, INTENT(INOUT) :: XAUX(NPOIN,*),AX(NELMAX,*)
00127       DOUBLE PRECISION, INTENT(INOUT) :: AD(NPOIN)
00128 !
00129       TYPE(BIEF_MESH) :: MESH
00130 !
00131 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00132 !
00133       INTEGER I1,I2,I3,I4,I5,I6,IELEM,IAN,ICOM,NPOIN2,IPLAN,K
00134       INTEGER S1,S2,S3,IT1,IT2,IELEM3D,L1,L2,ISEG
00135 !
00136 !     TETRA : WILL GIVE THE LOCAL NUMBERS OF POINTS IN THE PRISM
00137 !             THE 0 CORRESPOND TO SITUATIONS
00138 !             THAT NEVER HAPPEN (TETRA(1,1,1,... OR TETRA(2,2,2,...)
00139 !             SEE ALSO CPIKLE3 AND FLUX_EF_VF_3D, WITH SIMILAR USE
00140 !     INTEGER TETRA(2,2,2,3,4)
00141 !     DATA TETRA / 0,1,1,1,1,1,1,0,0,4,4,4,4,4,4,0,0,6,4,5,5,4,6,0,
00142 !    &             0,2,2,2,2,2,2,0,0,6,6,6,6,6,6,0,0,3,1,2,2,1,3,0,
00143 !    &             0,3,3,3,3,3,3,0,0,5,5,5,5,5,5,0,0,2,3,4,1,6,5,0,
00144 !    &             0,4,5,4,6,6,5,0,0,2,3,3,1,2,1,0,0,4,5,3,6,2,1,0 /
00145 !     EBE SYMMETRIC STORAGE OF OFF-DIAGONAL TERMS FOR TETRAHEDRONS
00146 !     0 ARE PUT FOR DIAGONALS
00147       INTEGER STO(4,4)
00148       DATA STO / 0, 1, 2, 3,
00149      &           1, 0, 4, 5,
00150      &           2, 4, 0, 6,
00151      &           3, 5, 6, 0 /
00152 !     CORRESPONDING OFF-DIAGONAL TERMS
00153 !     DATA STO /1-1,2-1,3-1,4-1,
00154 !    &          1-2,2-2,3-2,4-2,
00155 !    &          1-3,2-3,3-3,4-3,
00156 !    &          1-4,2-4,3-4,4-4/
00157 !     INTEGER ISEGT(6,2)
00158 !     DATA ISEGT/1,2,3,1,2,3,2,3,1,4,4,4/
00159 !
00160 !
00161 !-----------------------------------------------------------------------
00162 !
00163 !     CONSIDERS HERE THAT NPOIN < NELMAX TO USE XAUX AS XAUX(NPOIN,3)
00164 !
00165 !     XAUX(I,1) IS COEFFICIENT OF POINT BELOW I IN EQUATION OF POINT I
00166 !     XAUX(I,2) IS THE DIAGONAL
00167 !     XAUX(I,3) IS COEFFICIENT OF POINT ABOVE I IN EQUATION OF POINT I
00168 !
00169 !-----------------------------------------------------------------------
00170 !     INITIALISES THE DIAGONAL AND OFF-DIAGONAL TERMS
00171 !-----------------------------------------------------------------------
00172 !
00173       CALL OV('X=C     ',XAUX(1,1),AD,AD,0.D0,NPOIN)
00174       CALL OV('X=CY    ',XAUX(1,2),AD,AD,TETA,NPOIN)
00175       CALL OV('X=C     ',XAUX(1,3),AD,AD,0.D0,NPOIN)
00176 !
00177       CALL OV('X=CX    ',AD,AD,AD,1.D0-TETA,NPOIN)
00178 !
00179 !-----------------------------------------------------------------------
00180 !     ADDS TRIDIAGONAL TERMS
00181 !-----------------------------------------------------------------------
00182 !
00183       IF(IELM3.EQ.41) THEN
00184 !
00185         DO IELEM=1,NELEM
00186 !
00187           I1=IKLE(IELEM,1)
00188           I2=IKLE(IELEM,2)
00189           I3=IKLE(IELEM,3)
00190           I4=IKLE(IELEM,4)
00191           I5=IKLE(IELEM,5)
00192           I6=IKLE(IELEM,6)
00193           XAUX(I1,3)=XAUX(I1,3)+TETA*AX(IELEM,03) ! TERM 1-4
00194           XAUX(I2,3)=XAUX(I2,3)+TETA*AX(IELEM,08) ! TERM 2-5
00195           XAUX(I3,3)=XAUX(I3,3)+TETA*AX(IELEM,12) ! TERM 3-6
00196           XAUX(I4,1)=XAUX(I4,1)+TETA*AX(IELEM,03) ! TERM 4-1
00197           XAUX(I5,1)=XAUX(I5,1)+TETA*AX(IELEM,08) ! TERM 5-2
00198           XAUX(I6,1)=XAUX(I6,1)+TETA*AX(IELEM,12) ! TERM 6-3
00199 !
00200           AX(IELEM,03)=AX(IELEM,03)*(1.D0-TETA)
00201           AX(IELEM,08)=AX(IELEM,08)*(1.D0-TETA)
00202           AX(IELEM,12)=AX(IELEM,12)*(1.D0-TETA)
00203 !
00204         ENDDO
00205 !
00206       ELSEIF(IELM3.EQ.51) THEN
00207 !
00208         DO IPLAN=1,NPLAN-1
00209           DO IELEM=1,NELEM2
00210 !
00211 !           HERE LOWER LEVEL OF ELEMENTS ALWAYS TAKEN
00212 !           TO FIND THE WAY THE PRISM HAS BEEN CUT BY LOOKING
00213 !           AT GLOBAL NUMBERS OF POINTS
00214 !           THIS PART IS THUS COMMON TO ALL PLANES
00215 !           IKLE 3D IS TAKEN, COULD BE IKLE 2D AS WELL
00216             I1=IKLE(IELEM,1)
00217             I2=IKLE(IELEM,2)
00218             I3=IKLE(IELEM,3)
00219             IF(NCSIZE.GT.1) THEN
00220               I1=KNOLG(I1)
00221               I2=KNOLG(I2)
00222               I3=KNOLG(I3)
00223             ENDIF
00224 !           THIS IS DONE LIKE IN CPIKLE3 TO USE ARRAY TETRA
00225             IF(I1.GT.I2) THEN
00226               S1=1
00227             ELSE
00228               S1=2
00229             ENDIF
00230             IF(I2.GT.I3) THEN
00231               S2=1
00232             ELSE
00233               S2=2
00234             ENDIF
00235             IF(I3.GT.I1) THEN
00236               S3=1
00237             ELSE
00238               S3=2
00239             ENDIF
00240 !
00241 !           NOW TAKING CONTRIBUTIONS OF TETRAHEDRON K= 1, 2 AND 3
00242 !
00243             DO K=1,3
00244               IELEM3D=3*(IPLAN-1)*NELEM2+IELEM+(K-1)*NELEM2
00245 !             SEGMENTS 1 TO 6
00246               DO ISEG=1,6
00247 !               LOCAL NUMBERS OF 2 POINTS OF SEGMENT IN THE TETRAHEDRON
00248                 L1=ISEGT(ISEG,1)
00249                 L2=ISEGT(ISEG,2)
00250 !               GLOBAL NUMBERS OF 2 POINTS OF SEGMENT
00251                 I1=IKLE(IELEM3D,L1)
00252                 I2=IKLE(IELEM3D,L2)
00253 !               NUMBERS OF 2 POINTS OF SEGMENT IN THE ORIGINAL PRISM
00254                 IT1=TETRA(S1,S2,S3,K,L1)
00255                 IT2=TETRA(S1,S2,S3,K,L2)
00256 !               WE LOOK FOR VERTICALS OF THE ORIGINAL PRISM
00257                 IF(IT1.EQ.1.AND.IT2.EQ.4) THEN
00258                   XAUX(I1,3)=XAUX(I1,3)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 1-4
00259                   XAUX(I2,1)=XAUX(I2,1)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 4-1
00260                   AX(IELEM3D,STO(L1,L2))=
00261      &            AX(IELEM3D,STO(L1,L2))*(1.D0-TETA)
00262                 ELSEIF(IT1.EQ.4.AND.IT2.EQ.1) THEN
00263                   XAUX(I1,1)=XAUX(I1,1)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 4-1
00264                   XAUX(I2,3)=XAUX(I2,3)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 1-4
00265                   AX(IELEM3D,STO(L1,L2))=
00266      &            AX(IELEM3D,STO(L1,L2))*(1.D0-TETA)
00267                 ELSEIF(IT1.EQ.2.AND.IT2.EQ.5) THEN
00268                   XAUX(I1,3)=XAUX(I1,3)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 2-5
00269                   XAUX(I2,1)=XAUX(I2,1)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 5-2
00270                   AX(IELEM3D,STO(L1,L2))=
00271      &            AX(IELEM3D,STO(L1,L2))*(1.D0-TETA)
00272                 ELSEIF(IT1.EQ.5.AND.IT2.EQ.2) THEN
00273                   XAUX(I1,1)=XAUX(I1,1)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 5-2
00274                   XAUX(I2,3)=XAUX(I2,3)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 2-5
00275                   AX(IELEM3D,STO(L1,L2))=
00276      &            AX(IELEM3D,STO(L1,L2))*(1.D0-TETA)
00277                 ELSEIF(IT1.EQ.3.AND.IT2.EQ.6) THEN
00278                   XAUX(I1,3)=XAUX(I1,3)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 3-6
00279                   XAUX(I2,1)=XAUX(I2,1)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 6-3
00280                   AX(IELEM3D,STO(L1,L2))=
00281      &            AX(IELEM3D,STO(L1,L2))*(1.D0-TETA)
00282                 ELSEIF(IT1.EQ.6.AND.IT2.EQ.3) THEN
00283                   XAUX(I1,1)=XAUX(I1,1)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 6-3
00284                   XAUX(I2,3)=XAUX(I2,3)+TETA*AX(IELEM3D,STO(L1,L2)) ! TERM 3-6
00285                   AX(IELEM3D,STO(L1,L2))=
00286      &            AX(IELEM3D,STO(L1,L2))*(1.D0-TETA)
00287                 ENDIF
00288               ENDDO
00289             ENDDO
00290 !
00291         ENDDO
00292       ENDDO
00293 !
00294 !-----------------------------------------------------------------------
00295 !
00296       ELSE
00297         WRITE(LU,*) 'GETTRIEBE: ELEMENT NOT IMPLEMENTED:',IELM3
00298         CALL PLANTE(1)
00299         STOP
00300       ENDIF
00301 !
00302 !-----------------------------------------------------------------------
00303 !
00304 !     PARALLEL MODE
00305 !
00306       IF(NCSIZE.GT.1) THEN
00307         IAN    = 3
00308         ICOM   = 2
00309         NPOIN2 = BIEF_NBPTS(11,MESH)
00310         CALL PARCOM2(XAUX(1,1),XAUX(1,2),XAUX(1,3),
00311      &               NPOIN2,NPLAN,ICOM,IAN,MESH)
00312       ENDIF
00313 !
00314 !-----------------------------------------------------------------------
00315 !
00316       RETURN
00317       END

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