The TELEMAC-MASCARET system  trunk
mt06tt.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt06tt
3 ! *****************
4 !
5  &( t,xm,xmul,sf,f,x,y,z,ikle,nelem,nelmax)
6 !
7 !***********************************************************************
8 ! BIEF V6P1 21/08/2010
9 !***********************************************************************
10 !
11 !brief COMPUTES THE FMATMA MATRIX FOR TETRAHEDRONS.
12 !+
13 !+ THE VECTOR F CAN BE P0 OR P1.
14 !code
15 !+ STORAGE CONVENTION FOR EXTRA-DIAGONAL TERMS:
16 !+
17 !+ XM(IELEM, 1) ----> M(1,2) = M(2,1)
18 !+ XM(IELEM, 2) ----> M(1,3) = M(3,1)
19 !+ XM(IELEM, 3) ----> M(1,4) = M(4,1)
20 !+ XM(IELEM, 4) ----> M(2,3) = M(3,2)
21 !+ XM(IELEM, 5) ----> M(2,4) = M(4,2)
22 !+ XM(IELEM, 6) ----> M(3,4) = M(4,3)
23 !
24 !history J-M HERVOUET (LNH) ; F LEPEINTRE (LNH)
25 !+ 10/01/06
26 !+ V5P6
27 !+
28 !
29 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
30 !+ 13/07/2010
31 !+ V6P0
32 !+ Translation of French comments within the FORTRAN sources into
33 !+ English comments
34 !
35 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
36 !+ 21/08/2010
37 !+ V6P0
38 !+ Creation of DOXYGEN tags for automated documentation and
39 !+ cross-referencing of the FORTRAN sources
40 !
41 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
42 !| F |-->| FUNCTION USED IN THE FORMULA
43 !| FORMUL |-->| FORMULA DESCRIBING THE RESULTING MATRIX
44 !| IKLE |-->| CONNECTIVITY TABLE.
45 !| NELEM |-->| NUMBER OF ELEMENTS
46 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
47 !| SF |-->| STRUCTURE OF FUNCTIONS F
48 !| SURFAC |-->| AREA OF 2D ELEMENTS
49 !| T |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
50 !| Z |-->| ELEVATIONS OF POINTS
51 !| XM |<->| OFF-DIAGONAL TERMS
52 !| XMUL |-->| COEFFICIENT FOR MULTIPLICATION
53 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
54 !
55  USE bief, ex_mt06tt => mt06tt
57  IMPLICIT NONE
58 !
59 !
60 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
61 !
62  INTEGER, INTENT(IN) :: NELEM,NELMAX
63  INTEGER, INTENT(IN) :: IKLE(nelmax,4)
64 !
65  DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,4),XM(nelmax,6)
66 !
67  DOUBLE PRECISION, INTENT(IN) :: XMUL
68 !
69  DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*),Z(*)
70 !
71  TYPE(bief_obj), INTENT(IN) :: SF
72  DOUBLE PRECISION, INTENT(IN) :: F(*)
73 !
74 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
75 !
76 ! SPECIFIC DECLARATIONS
77 !
78  DOUBLE PRECISION X2,Y2,Z2,X3,Y3,Z3,X4,Y4,Z4,F1,F2,F3,F4
79  INTEGER I1,I2,I3,I4
80  INTEGER IELMF,IELEM
81 !
82  DOUBLE PRECISION XSUR720,JACOB
83 !
84 !***********************************************************************
85 !
86 ! DISCRETISES VECTOR F
87  ielmf=sf%ELM
88  xsur720=xmul/720.d0
89 !
90 !-----------------------------------------------------------------------
91 !
92 ! CASE WHERE F IS P1 LINEAR
93 !
94  IF (ielmf .EQ. 31) THEN
95 ! LOOP ON THE TETRAHEDRONS
96 !
97  DO ielem=1,nelem
98 !
99  i1=ikle(ielem,1)
100  i2=ikle(ielem,2)
101  i3=ikle(ielem,3)
102  i4=ikle(ielem,4)
103 !
104  f1 = f(i1)
105  f2 = f(i2)
106  f3 = f(i3)
107  f4 = f(i4)
108 !
109 !-----------------------------------------------------------------------
110 !
111  x2=x(i2)-x(i1)
112  y2=y(i2)-y(i1)
113  z2=z(i2)-z(i1)
114  x3=x(i3)-x(i1)
115  y3=y(i3)-y(i1)
116  z3=z(i3)-z(i1)
117  x4=x(i4)-x(i1)
118  y4=y(i4)-y(i1)
119  z4=z(i4)-z(i1)
120 !
121 ! JACOBIAN : Z2*(X3*Y4-X4*Y3)+Y2*(X4*Z3-X3*Z4)+X2*(Y3*Z4-Y4*Z3)
122 !
123 ! VOLUME OF THE TETRAHEDRON:
124 !
125 ! (Z2*(X3*Y4-X4*Y3)+Y2*(X4*Z3-X3*Z4)+X2*(Y3*Z4-Y4*Z3))/6
126 !
127  jacob = ( z2*(x3*y4-x4*y3)
128  & + y2*(x4*z3-x3*z4)
129  & + x2*(y3*z4-y4*z3) ) * xsur720
130 !
131  t(ielem,1) = jacob*(2*f3+6*f1+2*f4+2*f2)
132  t(ielem,2) = jacob*(2*f3+2*f1+2*f4+6*f2)
133  t(ielem,3) = jacob*2*(f1+3*f3+f2+f4)
134  t(ielem,4) = jacob*(2*f1+2*f3+6*f4+2*f2)
135 !
136  xm(ielem,1) = jacob*(f3+2*f1+f4+2*f2)
137  xm(ielem,2) = jacob*(2*f3+f2+2*f1+f4)
138  xm(ielem,3) = jacob*(f2+f3+2*f4+2*f1)
139  xm(ielem,4) = jacob*(f1+2*f2+f4+2*f3)
140  xm(ielem,5) = jacob*(2*f4+f3+f1+2*f2)
141  xm(ielem,6) = jacob*(f1+2*f3+f2+2*f4)
142 !
143 !-----------------------------------------------------------------------
144 !
145  ENDDO ! IELEM
146 !
147 !-----------------------------------------------------------------------
148 !
149 ! CASE WHERE F IS CONSTANT BY P0 ELEMENT
150 !
151  ELSEIF (ielmf .EQ. 30) THEN
152 ! LOOP ON THE TETRAHEDRONS
153 !
154  DO ielem=1,nelem
155 !
156  i1=ikle(ielem,1)
157  i2=ikle(ielem,2)
158  i3=ikle(ielem,3)
159  i4=ikle(ielem,4)
160 !
161 ! SAME VALUE FOR THE 4 FI
162 ! NOTE THAT THE COMPUTATIONS FOR T AND XM
163 ! COULD BE SIMPLIFIED ...
164 !
165  f1 = f(ielem)
166  f2 = f1
167  f3 = f1
168  f4 = f1
169 !
170 !-----------------------------------------------------------------------
171 !
172  x2=x(i2)-x(i1)
173  y2=y(i2)-y(i1)
174  z2=z(i2)-z(i1)
175  x3=x(i3)-x(i1)
176  y3=y(i3)-y(i1)
177  z3=z(i3)-z(i1)
178  x4=x(i4)-x(i1)
179  y4=y(i4)-y(i1)
180  z4=z(i4)-z(i1)
181 !
182 ! JACOBIAN : Z2*(X3*Y4-X4*Y3)+Y2*(X4*Z3-X3*Z4)+X2*(Y3*Z4-Y4*Z3)
183 !
184 ! VOLUME OF THE TETRAHEDRON:
185 !
186 ! (Z2*(X3*Y4-X4*Y3)+Y2*(X4*Z3-X3*Z4)+X2*(Y3*Z4-Y4*Z3))/6
187 !
188  jacob = ( z2*(x3*y4-x4*y3)
189  & + y2*(x4*z3-x3*z4)
190  & + x2*(y3*z4-y4*z3) ) * xsur720
191 !
192  t(ielem,1) = jacob*(2*f3+6*f1+2*f4+2*f2)
193  t(ielem,2) = jacob*(2*f3+2*f1+2*f4+6*f2)
194  t(ielem,3) = jacob*2*(f1+3*f3+f2+f4)
195  t(ielem,4) = jacob*(2*f1+2*f3+6*f4+2*f2)
196 !
197  xm(ielem,1) = jacob*(f3+2*f1+f4+2*f2)
198  xm(ielem,2) = jacob*(2*f3+f2+2*f1+f4)
199  xm(ielem,3) = jacob*(f2+f3+2*f4+2*f1)
200  xm(ielem,4) = jacob*(f1+2*f2+f4+2*f3)
201  xm(ielem,5) = jacob*(2*f4+f3+f1+2*f2)
202  xm(ielem,6) = jacob*(f1+2*f3+f2+2*f4)
203 !
204 !-----------------------------------------------------------------------
205 !
206  ENDDO ! IELEM
207 !
208 !-----------------------------------------------------------------------
209  ELSE
210 !
211  WRITE(lu,101) ielmf,sf%NAME
212 101 FORMAT(1x,'MT06TT (BIEF) :',/,
213  & 1x,'DISCRETIZATION OF F NOT AVAILABLE:',1i6,
214  & 1x,'REAL NAME: ',a6)
215  CALL plante(1)
216  stop
217 !
218  ENDIF
219 !
220 !-----------------------------------------------------------------------
221 !
222  RETURN
223  END SUBROUTINE mt06tt
subroutine mt06tt(T, XM, XMUL, SF, F, X, Y, Z, IKLE, NELEM, NELMAX)
Definition: mt06tt.f:7
Definition: bief.f:3