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