The TELEMAC-MASCARET system  trunk
mt08bb.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt08bb
3 ! *****************
4 !
5  &( a11 , a12 , a13 , a14 ,
6  & a21 , a22 , a23 , a24 ,
7  & a31 , a32 , a33 , a34 ,
8  & a41 , a42 , a43 , a44 ,
9  & xmul,sf,f,xel,yel,ikle1,ikle2,ikle3,ikle4,nelem,nelmax,icoord)
10 !
11 !***********************************************************************
12 ! BIEF V6P1 21/08/2010
13 !***********************************************************************
14 !
15 !brief COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
16 !code
17 !+ EXAMPLE WITH ICOORD = 1
18 !+
19 !+ / D
20 !+ A(I,J)=-XMUL / PSI2(J) * F * --( PSI1(I) ) D(OMEGA)
21 !+ /OMEGA DX
22 !+
23 !+ BEWARE THE MINUS SIGN !!
24 !+
25 !+ PSI1: BASES OF TYPE P1 TRIANGLE
26 !+ PSI2: BASES OF TYPE IELM2
27 !+
28 !+ IT WOULD BE A DERIVATIVE WRT Y WITH ICOORD=2
29 !
30 !warning THE JACOBIAN MUST BE POSITIVE
31 !
32 !history J-M HERVOUET (LNH) ; C MOULIN (LNH)
33 !+ 09/12/94
34 !+ V5P1
35 !+
36 !
37 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
38 !+ 13/07/2010
39 !+ V6P0
40 !+ Translation of French comments within the FORTRAN sources into
41 !+ English comments
42 !
43 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
44 !+ 21/08/2010
45 !+ V6P0
46 !+ Creation of DOXYGEN tags for automated documentation and
47 !+ cross-referencing of the FORTRAN sources
48 !
49 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
50 !| A11 |<--| ELEMENTS OF MATRIX
51 !| A12 |<--| ELEMENTS OF MATRIX
52 !| A13 |<--| ELEMENTS OF MATRIX
53 !| A14 |<--| ELEMENTS OF MATRIX
54 !| A21 |<--| ELEMENTS OF MATRIX
55 !| A22 |<--| ELEMENTS OF MATRIX
56 !| A23 |<--| ELEMENTS OF MATRIX
57 !| A24 |<--| ELEMENTS OF MATRIX
58 !| A31 |<--| ELEMENTS OF MATRIX
59 !| A32 |<--| ELEMENTS OF MATRIX
60 !| A33 |<--| ELEMENTS OF MATRIX
61 !| A34 |<--| ELEMENTS OF MATRIX
62 !| A41 |<--| ELEMENTS OF MATRIX
63 !| A42 |<--| ELEMENTS OF MATRIX
64 !| A43 |<--| ELEMENTS OF MATRIX
65 !| A44 |<--| ELEMENTS OF MATRIX
66 !| F |-->| FUNCTION USED IN THE FORMULA
67 !| ICOORD |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
68 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
69 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
70 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
71 !| NELEM |-->| NUMBER OF ELEMENTS
72 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
73 !| SF |-->| STRUCTURE OF FUNCTIONS F
74 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
75 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
76 !| XMUL |-->| MULTIPLICATION FACTOR
77 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
78 !
79  USE bief, ex_mt08bb => mt08bb
80 !
82  IMPLICIT NONE
83 !
84 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
85 !
86  INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
87  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax)
88  INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
89 !
90  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*),A14(*)
91  DOUBLE PRECISION, INTENT(INOUT) :: A21(*),A22(*),A23(*),A24(*)
92  DOUBLE PRECISION, INTENT(INOUT) :: A31(*),A32(*),A33(*),A34(*)
93  DOUBLE PRECISION, INTENT(INOUT) :: A41(*),A42(*),A43(*),A44(*)
94 !
95  DOUBLE PRECISION, INTENT(IN) :: XMUL
96  DOUBLE PRECISION, INTENT(IN) :: F(*)
97 !
98 ! STRUCTURE OF F
99  TYPE(bief_obj), INTENT(IN) :: SF
100 !
101  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
102 !
103 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
104 !
105  INTEGER IELEM,IELMF
106  DOUBLE PRECISION X2,X3,Y2,Y3,F1,F2,F3,F4
107 !
108 !-----------------------------------------------------------------------
109 !
110  ielmf=sf%ELM
111 !
112 !-----------------------------------------------------------------------
113 ! CASE WHERE F IS OF P1 DISCRETISATION
114 !-----------------------------------------------------------------------
115 !
116  IF(ielmf.EQ.11) THEN
117 !
118 !================================
119 ! CASE OF DERIVATIVE WRT X =
120 !================================
121 !
122  IF(icoord.EQ.1) THEN
123 !
124 ! LOOP ON THE ELEMENTS
125 !
126  DO ielem = 1 , nelem
127 !
128 ! INITIALISES THE GEOMETRICAL VARIABLES
129 !
130  y2 = yel(ielem,2)
131  y3 = yel(ielem,3)
132 !
133  f1 = f(ikle1(ielem))
134  f2 = f(ikle2(ielem))
135  f3 = f(ikle3(ielem))
136 !
137 ! EXTRADIAGONAL TERMS
138 !
139  a12(ielem) = (2*y2*(-f3-7*f2-4*f1)+y3*(f3+7*f2+4*f1))*xmul/216
140  a13(ielem) = (y2*(-7*f3-f2-4*f1)+2*y3*(7*f3+f2+4*f1))*xmul/216
141  a14(ielem) = (y2*(-3*f3-4*f2-5*f1)+y3*(4*f3+3*f2+5*f1))*xmul/72
142  a21(ielem) = (y2*(-f3-4*f2-7*f1)+y3*(-f3-4*f2-7*f1))*xmul/216
143  a23(ielem) = (y2*(7*f3+4*f2+f1)+2*y3*(-7*f3-4*f2-f1))*xmul/216
144  a24(ielem) = (y2*(f3-f1)+y3*(-4*f3-5*f2-3*f1))*xmul/72
145  a31(ielem) = (y2*(4*f3+f2+7*f1)+y3*(4*f3+f2+7*f1))*xmul/216
146  a32(ielem) = (2*y2*(4*f3+7*f2+f1)+y3*(-4*f3-7*f2-f1))*xmul/216
147  a34(ielem) = (y2*(5*f3+4*f2+3*f1)+y3*(-f2+f1))*xmul/72
148  a41(ielem) = (y2*(f3+4*f2+7*f1)+y3*(-4*f3-f2-7*f1))*xmul/72
149  a42(ielem) = (3*y2*(-f3+f1)+y3*(4*f3+7*f2+f1))*xmul/72
150  a43(ielem) = (y2*(-7*f3-4*f2-f1)+3*y3*(f2-f1))*xmul/72
151 !
152 ! DIAGONAL TERMS
153 ! (SUM OF EACH LINE IN THE MATRIX IS 0)
154 !
155  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
156  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
157  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
158  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
159 !
160  ENDDO ! IELEM
161 !
162  ELSEIF(icoord.EQ.2) THEN
163 !
164 !================================
165 ! CASE OF DERIVATIVE WRT Y =
166 !================================
167 !
168  DO ielem = 1 , nelem
169 !
170 ! INITIALISES THE GEOMETRICAL VARIABLES
171 !
172  x2 = xel(ielem,2)
173  x3 = xel(ielem,3)
174 !
175  f1 = f(ikle1(ielem))
176  f2 = f(ikle2(ielem))
177  f3 = f(ikle3(ielem))
178 !
179 ! EXTRADIAGONAL TERMS
180 !
181  a12(ielem) = (2*x2*(f3+7*f2+4*f1)+x3*(-f3-7*f2-4*f1))*xmul/216
182  a13(ielem) = (x2*(7*f3+f2+4*f1)+2*x3*(-7*f3-f2-4*f1))*xmul/216
183  a14(ielem) = (x2*(3*f3+4*f2+5*f1)+x3*(-4*f3-3*f2-5*f1))*xmul/72
184  a21(ielem) = (x2*(f3+4*f2+7*f1)+x3*(f3+4*f2+7*f1))*xmul/216
185  a23(ielem) = (x2*(-7*f3-4*f2-f1)+2*x3*(7*f3+4*f2+f1))*xmul/216
186  a24(ielem) = (x2*(-f3+f1)+x3*(4*f3+5*f2+3*f1))*xmul/72
187  a31(ielem) = (x2*(-4*f3-f2-7*f1)+x3*(-4*f3-f2-7*f1))*xmul/216
188  a32(ielem) = (2*x2*(-4*f3-7*f2-f1)+x3*(4*f3+7*f2+f1))*xmul/216
189  a34(ielem) = (x2*(-5*f3-4*f2-3*f1)+x3*(f2-f1))*xmul/72
190  a41(ielem) = (x2*(-f3-4*f2-7*f1)+x3*(4*f3+f2+7*f1))*xmul/72
191  a42(ielem) = (3*x2*(f3-f1)+x3*(-4*f3-7*f2-f1))*xmul/72
192  a43(ielem) = (x2*(7*f3+4*f2+f1)+3*x3*(-f2+f1))*xmul/72
193 !
194 ! DIAGONAL TERMS
195 ! (SUM OF EACH LINE IN THE MATRIX IS 0)
196 !
197  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
198  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
199  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
200  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
201 !
202  ENDDO ! IELEM
203 !
204  ELSE
205 !
206  WRITE(lu,201) icoord
207  CALL plante(0)
208  stop
209  ENDIF
210 !
211 !
212 !-----------------------------------------------------------------------
213 ! CASE WHERE F IS OF QUASI-BUBBLE DISCRETISATION
214 !-----------------------------------------------------------------------
215 !
216  ELSEIF(ielmf.EQ.12) THEN
217 !
218 !================================
219 ! CASE OF DERIVATIVE WRT X =
220 !================================
221 !
222  IF(icoord.EQ.1) THEN
223 !
224 ! LOOP ON THE ELEMENTS
225 !
226  DO ielem = 1 , nelem
227 !
228 ! INITIALISES THE GEOMETRICAL VARIABLES
229 !
230  y2 = yel(ielem,2)
231  y3 = yel(ielem,3)
232 !
233  f1 = f(ikle1(ielem))
234  f2 = f(ikle2(ielem))
235  f3 = f(ikle3(ielem))
236  f4 = f(ikle4(ielem))
237 !
238 ! EXTRADIAGONAL TERMS
239 !
240  a12(ielem) = (2*y2*(-f4-2*f2-f1)+y3*(f4+2*f2+f1))*xmul/72
241  a13(ielem) = (y2*(-2*f3-f4-f1)+2*y3*(2*f3+f4+f1))*xmul/72
242  a14(ielem) = (y2*(-f3-6*f4-2*f2-3*f1)+y3*(2*f3+6*f4+f2+
243  & 3*f1))*xmul/72
244  a21(ielem) = (y2*(-f4-f2-2*f1)+y3*(-f4-f2-2*f1))*xmul/72
245  a23(ielem) = (y2*(2*f3+f4+f2)+2*y3*(-2*f3-f4-f2))*xmul/72
246  a24(ielem) = (y2*(f3-f1)+y3*(-2*f3-6*f4-3*f2-f1))*xmul/72
247  a31(ielem) = (y2*(f3+f4+2*f1)+y3*(f3+f4+2*f1))*xmul/72
248  a32(ielem) = (2*y2*(f3+f4+2*f2)+y3*(-f3-f4-2*f2))*xmul/72
249  a34(ielem) = (y2*(3*f3+6*f4+2*f2+f1)+y3*(-f2+f1))*xmul/72
250  a41(ielem) = (y2*(f4+f2+2*f1)+y3*(-f3-f4-2*f1))*xmul/24
251  a42(ielem) = (y2*(-f3+f1)+y3*(f3+f4+2*f2))*xmul/24
252  a43(ielem) = (y2*(-2*f3-f4-f2)+y3*(f2-f1))*xmul/24
253 !
254 ! DIAGONAL TERMS
255 ! (SUM OF EACH LINE IN THE MATRIX IS 0)
256 !
257  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
258  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
259  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
260  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
261 !
262  ENDDO ! IELEM
263 !
264  ELSEIF(icoord.EQ.2) THEN
265 !
266 !================================
267 ! CASE OF DERIVATIVE WRT Y =
268 !================================
269 !
270  DO ielem = 1 , nelem
271 !
272 ! INITIALISES THE GEOMETRICAL VARIABLES
273 !
274  x2 = xel(ielem,2)
275  x3 = xel(ielem,3)
276 !
277  f1 = f(ikle1(ielem))
278  f2 = f(ikle2(ielem))
279  f3 = f(ikle3(ielem))
280  f4 = f(ikle4(ielem))
281 !
282 ! EXTRADIAGONAL TERMS
283 !
284  a12(ielem) = (2*x2*(f4+2*f2+f1)+x3*(-f4-2*f2-f1))*xmul/72
285  a13(ielem) = (x2*(2*f3+f4+f1)+2*x3*(-2*f3-f4-f1))*xmul/72
286  a14(ielem) = (x2*(f3+6*f4+2*f2+3*f1)+x3*(-2*f3-6*f4-f2-
287  & 3*f1))*xmul/72
288  a21(ielem) = (x2*(f4+f2+2*f1)+x3*(f4+f2+2*f1))*xmul/72
289  a23(ielem) = (x2*(-2*f3-f4-f2)+2*x3*(2*f3+f4+f2))*xmul/72
290  a24(ielem) = (x2*(-f3+f1)+x3*(2*f3+6*f4+3*f2+f1))*xmul/72
291  a31(ielem) = (x2*(-f3-f4-2*f1)+x3*(-f3-f4-2*f1))*xmul/72
292  a32(ielem) = (2*x2*(-f3-f4-2*f2)+x3*(f3+f4+2*f2))*xmul/72
293  a34(ielem) = (x2*(-3*f3-6*f4-2*f2-f1)+x3*(f2-f1))*xmul/72
294  a41(ielem) = (x2*(-f4-f2-2*f1)+x3*(f3+f4+2*f1))*xmul/24
295  a42(ielem) = (x2*(f3-f1)+x3*(-f3-f4-2*f2))*xmul/24
296  a43(ielem) = (x2*(2*f3+f4+f2)+x3*(-f2+f1))*xmul/24
297 !
298 ! DIAGONAL TERMS
299 ! (SUM OF EACH LINE IN THE MATRIX IS 0)
300 !
301  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
302  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
303  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
304  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
305 !
306  ENDDO ! IELEM
307 !
308  ELSE
309 !
310  WRITE(lu,201) icoord
311  CALL plante(0)
312  stop
313  ENDIF
314 !
315 !-----------------------------------------------------------------------
316 !
317  ELSE
318  WRITE(lu,101) ielmf
319 101 FORMAT(1x,'MT08BB (BIEF) :',/,
320  & 1x,'DISCRETIZATION OF F : ',1i6,' NOT AVAILABLE')
321  CALL plante(0)
322  stop
323  ENDIF
324 !
325 201 FORMAT(1x,'MT08BB (BIEF) : IMPOSSIBLE COMPONENT ',
326  & 1i6,' CHECK ICOORD')
327 !
328 !-----------------------------------------------------------------------
329 !
330  RETURN
331  END
subroutine mt08bb(A11, A12, A13, A14, A21, A22, A23, A24, A31, A32, A33, A34, A41, A42, A43, A44, XMUL, SF, F, XEL, YEL, IKLE1, IKLE2, IKLE3, IKLE4, NELEM, NELMAX, ICOORD)
Definition: mt08bb.f:11
Definition: bief.f:3