The TELEMAC-MASCARET system  trunk
mt06pp.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt06pp
3 ! *****************
4 !
5  &( t,xm,xmul,sf,f,z,surfac,ikle,nelem,nelmax)
6 !
7 !***********************************************************************
8 ! BIEF V6P1 21/08/2010
9 !***********************************************************************
10 !
11 !brief BUILDS THE MASS MATRIX FOR P1 PRISMS.
12 !code
13 !+ COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
14 !+
15 !+ /
16 !+ A = / F * (P *P )*J(X,Y) DXDY
17 !+ I J /S I J
18 !+
19 !+ BY ELEMENTARY CELL; THE ELEMENT IS THE P1 PRISM
20 !+
21 !+ J(X,Y): JACOBIAN OF THE ISOPARAMETRIC TRANSFORMATION
22 !
23 !warning THE JACOBIAN MUST BE POSITIVE
24 !
25 !history J-M HERVOUET (LNH)
26 !+ 28/11/94
27 !+
28 !+
29 !
30 !history ARNAUD DESITTER - UNIVERSITY OF BRISTOL
31 !+ **/04/98
32 !+ V5P1
33 !+
34 !
35 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
36 !+ 13/07/2010
37 !+ V6P0
38 !+ Translation of French comments within the FORTRAN sources into
39 !+ English comments
40 !
41 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
42 !+ 21/08/2010
43 !+ V6P0
44 !+ Creation of DOXYGEN tags for automated documentation and
45 !+ cross-referencing of the FORTRAN sources
46 !
47 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
48 !| F |-->| FUNCTION USED IN THE FORMULA
49 !| FORMUL |-->| FORMULA DESCRIBING THE RESULTING MATRIX
50 !| IKLE |-->| CONNECTIVITY TABLE.
51 !| NELEM |-->| NUMBER OF ELEMENTS
52 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
53 !| SF |-->| STRUCTURE OF FUNCTIONS F
54 !| SURFAC |-->| AREA OF 2D ELEMENTS
55 !| T |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
56 !| Z |-->| ELEVATIONS OF POINTS
57 !| XM |<->| OFF-DIAGONAL TERMS
58 !| XMUL |-->| COEFFICIENT FOR MULTIPLICATION
59 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
60 !
61  USE bief, ex_mt06pp => mt06pp
62  USE declarations_special
63  IMPLICIT NONE
64 !
65 !
66 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
67 !
68  INTEGER, INTENT(IN) :: NELEM,NELMAX
69  INTEGER, INTENT(IN) :: IKLE(nelmax,6)
70 !
71  DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,6), XM(nelmax,30)
72  DOUBLE PRECISION, INTENT(IN) :: XMUL
73 !
74  DOUBLE PRECISION, INTENT(IN) :: F(*)
75 !
76 ! STRUCTURE OF F
77 !
78  TYPE(bief_obj), INTENT(IN) :: SF
79 !
80  DOUBLE PRECISION, INTENT(IN) :: Z(*)
81  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
82 !
83 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
84 !
85 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
86 !
87  DOUBLE PRECISION SUR2160,COEF,H1,H2,H3,F1,F2,F3,F4,F5,F6
88 !
89  INTEGER IELEM,IELMF
90 !
91 !-----------------------------------------------------------------------
92 !
93  sur2160 = xmul / 2160.d0
94 !
95 !-----------------------------------------------------------------------
96 !
97  ielmf=sf%ELM
98  IF(ielmf.NE.41) THEN
99  WRITE(lu,101) ielmf
100 101 FORMAT(1x,'MT06PP (BIEF) :',/,
101  & 1x,'DISCRETIZATION OF F : ',1i6,' NOT AVAILABLE')
102  CALL plante(1)
103  stop
104  ENDIF
105 !
106 !-----------------------------------------------------------------------
107 !
108 ! LOOP ON THE ELEMENTS
109 !
110  DO ielem = 1,nelem
111 !
112  coef = surfac(ielem)* sur2160
113 !
114  h1 = (z(ikle(ielem,4)) - z(ikle(ielem,1))) * coef
115  h2 = (z(ikle(ielem,5)) - z(ikle(ielem,2))) * coef
116  h3 = (z(ikle(ielem,6)) - z(ikle(ielem,3))) * coef
117 !
118  f1 = f(ikle(ielem,1))
119  f2 = f(ikle(ielem,2))
120  f3 = f(ikle(ielem,3))
121  f4 = f(ikle(ielem,4))
122  f5 = f(ikle(ielem,5))
123  f6 = f(ikle(ielem,6))
124 !
125 !-----------------------------------------------------------------------
126 !
127 ! EXTRA-DIAGONAL TERMS
128 !
129  xm(ielem,01) =
130  & (9*f1+6*f2+3*f3+3*f4+2*f5+f6)*h1+
131  & (6*f1+9*f2+3*f3+2*f4+3*f5+f6)*h2+
132  & (3*f1+3*f2+3*f3+f4+f5+f6)*h3
133  xm(ielem,02) =
134  & (9*f1+3*f2+6*f3+3*f4+f5+2*f6)*h1+
135  & (3*f1+3*f2+3*f3+f4+f5+f6)*h2+
136  & (6*f1+3*f2+9*f3+2*f4+f5+3*f6)*h3
137  xm(ielem,03) =
138  & (12*f1+3*f2+3*f3+12*f4+3*f5+3*f6)*h1+
139  & (3*f1+2*f2+f3+3*f4+2*f5+f6)*h2+
140  & (3*f1+f2+2*f3+3*f4+f5+2*f6)*h3
141  xm(ielem,04) =
142  & (3*f1+2*f2+f3+3*f4+2*f5+f6)*h1+
143  & (2*f1+3*f2+f3+2*f4+3*f5+f6)*h2+
144  & (f1+f2+f3+f4+f5+f6)*h3
145  xm(ielem,05) =
146  & (3*f1+f2+2*f3+3*f4+f5+2*f6)*h1+
147  & (f1+f2+f3+f4+f5+f6)*h2+
148  & (2*f1+f2+3*f3+2*f4+f5+3*f6)*h3
149  xm(ielem,06) =
150  & (3*f1+3*f2+3*f3+f4+f5+f6)*h1+
151  & (3*f1+9*f2+6*f3+f4+3*f5+2*f6)*h2+
152  & (3*f1+6*f2+9*f3+f4+2*f5+3*f6)*h3
153  xm(ielem,07) =
154  & (3*f1+2*f2+f3+3*f4+2*f5+f6)*h1+
155  & (2*f1+3*f2+f3+2*f4+3*f5+f6)*h2+
156  & (f1+f2+f3+f4+f5+f6)*h3
157  xm(ielem,08) =
158  & (2*f1+3*f2+f3+2*f4+3*f5+f6)*h1+
159  & (3*f1+12*f2+3*f3+3*f4+12*f5+3*f6)*h2+
160  & (f1+3*f2+2*f3+f4+3*f5+2*f6)*h3
161  xm(ielem,09) =
162  & (f1+f2+f3+f4+f5+f6)*h1+
163  & (f1+3*f2+2*f3+f4+3*f5+2*f6)*h2+
164  & (f1+2*f2+3*f3+f4+2*f5+3*f6)*h3
165  xm(ielem,10) =
166  & (3*f1+f2+2*f3+3*f4+f5+2*f6)*h1+
167  & (f1+f2+f3+f4+f5+f6)*h2+
168  & (2*f1+f2+3*f3+2*f4+f5+3*f6)*h3
169  xm(ielem,11) =
170  & (f1+f2+f3+f4+f5+f6)*h1+
171  & (f1+3*f2+2*f3+f4+3*f5+2*f6)*h2+
172  & (f1+2*f2+3*f3+f4+2*f5+3*f6)*h3
173  xm(ielem,12) =
174  & (2*f1+f2+3*f3+2*f4+f5+3*f6)*h1+
175  & (f1+2*f2+3*f3+f4+2*f5+3*f6)*h2+
176  & (3*f1+3*f2+12*f3+3*f4+3*f5+12*f6)*h3
177  xm(ielem,13) =
178  & (3*f1+2*f2+f3+9*f4+6*f5+3*f6)*h1+
179  & (2*f1+3*f2+f3+6*f4+9*f5+3*f6)*h2+
180  & (f1+f2+f3+3*f4+3*f5+3*f6)*h3
181  xm(ielem,14) =
182  & (3*f1+f2+2*f3+9*f4+3*f5+6*f6)*h1+
183  & (f1+f2+f3+3*f4+3*f5+3*f6)*h2+
184  & (2*f1+f2+3*f3+6*f4+3*f5+9*f6)*h3
185  xm(ielem,15) =
186  & (f1+f2+f3+3*f4+3*f5+3*f6)*h1+
187  & (f1+3*f2+2*f3+3*f4+9*f5+6*f6)*h2+
188  & (f1+2*f2+3*f3+3*f4+6*f5+9*f6)*h3
189 !
190 ! DIAGONAL TERMS
191 !
192  t(ielem,1) =
193  & (36*f1+9*f2+9*f3+12*f4+3*f5+3*f6)*h1+
194  & (9*f1+6*f2+3*f3+3*f4+2*f5+f6)*h2+
195  & (9*f1+3*f2+6*f3+3*f4+f5+2*f6)*h3
196  t(ielem,2) =
197  & (6*f1+9*f2+3*f3+2*f4+3*f5+f6)*h1+
198  & (9*f1+36*f2+9*f3+3*f4+12*f5+3*f6)*h2+
199  & (3*f1+9*f2+6*f3+f4+3*f5+2*f6)*h3
200  t(ielem,3) =
201  & (6*f1+3*f2+9*f3+2*f4+f5+3*f6)*h1+
202  & (3*f1+6*f2+9*f3+f4+2*f5+3*f6)*h2+
203  & (9*f1+9*f2+36*f3+3*f4+3*f5+12*f6)*h3
204  t(ielem,4) =
205  & (12*f1+3*f2+3*f3+36*f4+9*f5+9*f6)*h1+
206  & (3*f1+2*f2+f3+9*f4+6*f5+3*f6)*h2+
207  & (3*f1+f2+2*f3+9*f4+3*f5+6*f6)*h3
208  t(ielem,5) =
209  & (2*f1+3*f2+f3+6*f4+9*f5+3*f6)*h1+
210  & (3*f1+12*f2+3*f3+9*f4+36*f5+9*f6)*h2+
211  & (f1+3*f2+2*f3+3*f4+9*f5+6*f6)*h3
212  t(ielem,6) =
213  & (2*f1+f2+3*f3+6*f4+3*f5+9*f6)*h1+
214  & (f1+2*f2+3*f3+3*f4+6*f5+9*f6)*h2+
215  & (3*f1+3*f2+12*f3+9*f4+9*f5+36*f6)*h3
216 !
217  ENDDO
218 !
219 !-----------------------------------------------------------------------
220 !
221  RETURN
222  END
mt06pp
subroutine mt06pp(T, XM, XMUL, SF, F, Z, SURFAC, IKLE, NELEM, NELMAX)
Definition: mt06pp.f:7
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:09 for The TELEMAC-MASCARET system by doxygen 1.8.13