The TELEMAC-MASCARET system  trunk
mt03bb.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt03bb
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,sg,su,sv,f,g,u,v,
10  & xel,yel,ikle1,ikle2,ikle3,nelem,nelmax)
11 !
12 !***********************************************************************
13 ! BIEF V6P1 21/08/2010
14 !***********************************************************************
15 !
16 !brief COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
17 !code
18 !+ / --> - --> --> --->
19 !+ A(I,J) = XMUL / KEL . GRAD(PSI1(I)) * U . GRAD(PSI2(J)) D(OMEGA)
20 !+ /OMEGA
21 !+ -->
22 !+ KEL CONSTANT VECTOR ON THE ELEMENT, WITH COMPONENTS F AND G
23 !+
24 !+ PSI1 OF P1 DISCRETISATION
25 !+ PSI2 OF P2 DISCRETISATION
26 !
27 !warning THE JACOBIAN MUST BE POSITIVE
28 !
29 !history J-M HERVOUET (LNH) ; C MOULIN (LNH)
30 !+ 12/04/93
31 !+ V5P1
32 !+
33 !
34 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
35 !+ 13/07/2010
36 !+ V6P0
37 !+ Translation of French comments within the FORTRAN sources into
38 !+ English comments
39 !
40 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
41 !+ 21/08/2010
42 !+ V6P0
43 !+ Creation of DOXYGEN tags for automated documentation and
44 !+ cross-referencing of the FORTRAN sources
45 !
46 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
47 !| A11 |<--| ELEMENTS OF MATRIX
48 !| A12 |<--| ELEMENTS OF MATRIX
49 !| A13 |<--| ELEMENTS OF MATRIX
50 !| A14 |<--| ELEMENTS OF MATRIX
51 !| A21 |<--| ELEMENTS OF MATRIX
52 !| A22 |<--| ELEMENTS OF MATRIX
53 !| A23 |<--| ELEMENTS OF MATRIX
54 !| A24 |<--| ELEMENTS OF MATRIX
55 !| A31 |<--| ELEMENTS OF MATRIX
56 !| A32 |<--| ELEMENTS OF MATRIX
57 !| A33 |<--| ELEMENTS OF MATRIX
58 !| A34 |<--| ELEMENTS OF MATRIX
59 !| A41 |<--| ELEMENTS OF MATRIX
60 !| A42 |<--| ELEMENTS OF MATRIX
61 !| A43 |<--| ELEMENTS OF MATRIX
62 !| A44 |<--| ELEMENTS OF MATRIX
63 !| F |-->| FUNCTION USED IN THE FORMULA
64 !| G |-->| FUNCTION USED IN THE FORMULA
65 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
66 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
67 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
68 !| NELEM |-->| NUMBER OF ELEMENTS
69 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
70 !| SF |-->| STRUCTURE OF FUNCTIONS F
71 !| SG |-->| STRUCTURE OF FUNCTIONS G
72 !| SU |-->| BIEF_OBJ STRUCTURE OF U
73 !| SV |-->| BIEF_OBJ STRUCTURE OF V
74 !| U |-->| FUNCTION U USED IN THE FORMULA
75 !| V |-->| FUNCTION V USED IN THE FORMULA
76 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
77 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
78 !| XMUL |-->| MULTIPLICATION FACTOR
79 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
80 !
81  USE bief, ex_mt03bb => mt03bb
82 !
83  USE declarations_special
84  IMPLICIT NONE
85 !
86 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
87 !
88  INTEGER, INTENT(IN) :: NELEM,NELMAX
89  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax),IKLE3(nelmax)
90 !
91  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*),A14(*)
92  DOUBLE PRECISION, INTENT(INOUT) :: A21(*),A22(*),A23(*),A24(*)
93  DOUBLE PRECISION, INTENT(INOUT) :: A31(*),A32(*),A33(*),A34(*)
94  DOUBLE PRECISION, INTENT(INOUT) :: A41(*),A42(*),A43(*),A44(*)
95 !
96  DOUBLE PRECISION, INTENT(IN) :: XMUL
97  DOUBLE PRECISION, INTENT(IN) :: F(*),G(*),U(*),V(*)
98 !
99 ! STRUCTURES OF F, G, U, V
100  TYPE(bief_obj), INTENT(IN) :: SF,SG,SU,SV
101 !
102  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
103 !
104 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
105 !
106 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
107 !
108  INTEGER IELEM,IELMF,IELMG,IELMU,IELMV
109 !
110  DOUBLE PRECISION X2,X3,Y2,Y3,U1,U2,U3,U4,V1,V2,V3,V4
111  DOUBLE PRECISION KXEL,KYEL
112 !
113 !-----------------------------------------------------------------------
114 !
115  ielmf = sf%ELM
116  ielmg = sg%ELM
117  ielmu = su%ELM
118  ielmv = sv%ELM
119 !
120 !-----------------------------------------------------------------------
121 ! CASE WHERE U IS OF P1 DISCRETISATION
122 !-----------------------------------------------------------------------
123 !
124  IF(ielmf.EQ.10.AND.ielmg.EQ.10.AND.
125  & ielmu.EQ.11.AND.ielmv.EQ.11) THEN
126 !
127 ! LOOP ON THE ELEMENTS
128 !
129  DO ielem = 1 , nelem
130 !
131  x2 = xel(ielem,2)
132  x3 = xel(ielem,3)
133  y2 = yel(ielem,2)
134  y3 = yel(ielem,3)
135 !
136  u1 = u(ikle1(ielem))
137  u2 = u(ikle2(ielem))
138  u3 = u(ikle3(ielem))
139  v1 = v(ikle1(ielem))
140  v2 = v(ikle2(ielem))
141  v3 = v(ikle3(ielem))
142 !
143  kxel = f(ielem)
144  kyel = g(ielem)
145 !
146 ! COMPUTES 9 OF THE 16 TERMS
147 !
148  a11(ielem)=
149  & (x2**2*kyel*(8*v3+17*v2+20*v1)+4*x2*x3*kyel*(-
150  & 5*v3-5*v2-8*v1)+x2*y2*kxel*(-8*v3-17*v2-20*v1)+x2*y2*
151  & kyel*(-8*u3-17*u2-20*u1)+2*x2*y3*kxel*(5*v3+5*v2+8*v1)
152  & +2*x2*y3*kyel*(5*u3+5*u2+8*u1)+x3**2*kyel*(17*v3+8*v2+
153  & 20*v1)+2*x3*y2*kxel*(5*v3+5*v2+8*v1)+2*x3*y2*kyel*(5*u3
154  & +5*u2+8*u1)+x3*y3*kxel*(-17*v3-8*v2-20*v1)+x3*y3*kyel*(-
155  & 17*u3-8*u2-20*u1)+y2**2*kxel*(8*u3+17*u2+20*u1)+4*y2*
156  & y3*kxel*(-5*u3-5*u2-8*u1)+y3**2*kxel*(17*u3+8*u2+20*u1)
157  & )*xmul/(54*x2*y3-54*x3*y2)
158 !
159  a12(ielem)=
160  & (2*x2**2*kyel*(v3+4*v2+4*v1)+x2*x3*kyel*(v3+4*
161  & v2+4*v1)+2*x2*y2*kxel*(-v3-4*v2-4*v1)+2*x2*y2*kyel*(-u3-
162  & 4*u2-4*u1)+x2*y3*kxel*(v3+4*v2+4*v1)+2*x2*y3*kyel*(-u3-4
163  & *u2-4*u1)+x3**2*kyel*(-v3-4*v2-4*v1)+2*x3*y2*kxel*(-v3-4
164  & *v2-4*v1)+x3*y2*kyel*(u3+4*u2+4*u1)+x3*y3*kxel*(v3+4*v2+
165  & 4*v1)+x3*y3*kyel*(u3+4*u2+4*u1)+2*y2**2*kxel*(u3+4*u2+4*
166  & u1)+y2*y3*kxel*(u3+4*u2+4*u1)+y3**2*kxel*(-u3-4*u2-4*u1))
167  & *xmul/(54*x2*y3-54*x3*y2)
168 !
169  a13(ielem)=
170  & (x2**2*kyel*(-4*v3-v2-4*v1)+x2*x3*kyel*(4*v3+v2+
171  & 4*v1)+x2*y2*kxel*(4*v3+v2+4*v1)+x2*y2*kyel*(4*u3+u2+4*u1)
172  & +2*x2*y3*kxel*(-4*v3-v2-4*v1)+x2*y3*kyel*(4*u3+u2+4*u1)+
173  & 2*x3**2*kyel*(4*v3+v2+4*v1)+x3*y2*kxel*(4*v3+v2+4*v1)+2*
174  & x3*y2*kyel*(-4*u3-u2-4*u1)+2*x3*y3*kxel*(-4*v3-v2-4*v1)+
175  & 2*x3*y3*kyel*(-4*u3-u2-4*u1)+y2**2*kxel*(-4*u3-u2-4*u1)+
176  & y2*y3*kxel*(4*u3+u2+4*u1)+2*y3**2*kxel*(4*u3+u2+4*u1))
177  & *xmul/(54*x2*y3-54*x3*y2)
178 !
179  a21(ielem)=
180  & (2*x2**2*kyel*(v3+4*v2+4*v1)+x2*x3*kyel*(v3+4*
181  & v2+4*v1)+2*x2*y2*kxel*(-v3-4*v2-4*v1)+2*x2*y2*kyel*(-u3-
182  & 4*u2-4*u1)+2*x2*y3*kxel*(-v3-4*v2-4*v1)+x2*y3*kyel*(u3+4
183  & *u2+4*u1)+x3**2*kyel*(-v3-4*v2-4*v1)+x3*y2*kxel*(v3+4*v2+
184  & 4*v1)+2*x3*y2*kyel*(-u3-4*u2-4*u1)+x3*y3*kxel*(v3+4*v2+4
185  & *v1)+x3*y3*kyel*(u3+4*u2+4*u1)+2*y2**2*kxel*(u3+4*u2+4*
186  & u1)+y2*y3*kxel*(u3+4*u2+4*u1)+y3**2*kxel*(-u3-4*u2-4*u1))
187  & *xmul/(54*x2*y3-54*x3*y2)
188 !
189  a23(ielem)=
190  & (2*x2**2*kyel*(4*v3+4*v2+v1)+5*x2*x3*kyel*(-4*
191  & v3-4*v2-v1)+2*x2*y2*kxel*(-4*v3-4*v2-v1)+2*x2*y2*kyel*(-
192  & 4*u3-4*u2-u1)+4*x2*y3*kxel*(4*v3+4*v2+v1)+x2*y3*kyel*(4*
193  & u3+4*u2+u1)+2*x3**2*kyel*(4*v3+4*v2+v1)+x3*y2*kxel*(4*v3
194  & +4*v2+v1)+4*x3*y2*kyel*(4*u3+4*u2+u1)+2*x3*y3*kxel*(-4*
195  & v3-4*v2-v1)+2*x3*y3*kyel*(-4*u3-4*u2-u1)+2*y2**2*kxel*(
196  & 4*u3+4*u2+u1)+5*y2*y3*kxel*(-4*u3-4*u2-u1)+2*y3**2*kxel*
197  & (4*u3+4*u2+u1))*xmul/(54*x2*y3-54*x3*y2)
198 !
199  a24(ielem)=
200  & (x2**2*kyel*(-5*v3-8*v2-5*v1)+x2*x3*kyel*(11*v3
201  & +8*v2-v1)+x2*y2*kxel*(5*v3+8*v2+5*v1)+x2*y2*kyel*(5*u3+
202  & 8*u2+5*u1)+x2*y3*kxel*(-7*v3-4*v2+2*v1)+x2*y3*kyel*(-4*
203  & u3-4*u2-u1)+2*x3**2*kyel*(-4*v3-4*v2-v1)+x3*y2*kxel*(-4*
204  & v3-4*v2-v1)+x3*y2*kyel*(-7*u3-4*u2+2*u1)+2*x3*y3*kxel*(
205  & 4*v3+4*v2+v1)+2*x3*y3*kyel*(4*u3+4*u2+u1)+y2**2*kxel*(-5
206  & *u3-8*u2-5*u1)+y2*y3*kxel*(11*u3+8*u2-u1)+2*y3**2*kxel*(
207  & -4*u3-4*u2-u1))*xmul/(18*x2*y3-18*x3*y2)
208 !
209  a31(ielem)=
210  & (x2**2*kyel*(-4*v3-v2-4*v1)+x2*x3*kyel*(4*v3+v2+
211  & 4*v1)+x2*y2*kxel*(4*v3+v2+4*v1)+x2*y2*kyel*(4*u3+u2+4*u1)
212  & +x2*y3*kxel*(4*v3+v2+4*v1)+2*x2*y3*kyel*(-4*u3-u2-4*u1)+
213  & 2*x3**2*kyel*(4*v3+v2+4*v1)+2*x3*y2*kxel*(-4*v3-v2-4*v1)
214  & +x3*y2*kyel*(4*u3+u2+4*u1)+2*x3*y3*kxel*(-4*v3-v2-4*v1)+
215  & 2*x3*y3*kyel*(-4*u3-u2-4*u1)+y2**2*kxel*(-4*u3-u2-4*u1)+
216  & y2*y3*kxel*(4*u3+u2+4*u1)+2*y3**2*kxel*(4*u3+u2+4*u1))
217  & *xmul/(54*x2*y3-54*x3*y2)
218 !
219  a32(ielem)=
220  & (2*x2**2*kyel*(4*v3+4*v2+v1)+5*x2*x3*kyel*(-4*
221  & v3-4*v2-v1)+2*x2*y2*kxel*(-4*v3-4*v2-v1)+2*x2*y2*kyel*(-
222  & 4*u3-4*u2-u1)+x2*y3*kxel*(4*v3+4*v2+v1)+4*x2*y3*kyel*(4*
223  & u3+4*u2+u1)+2*x3**2*kyel*(4*v3+4*v2+v1)+4*x3*y2*kxel*(4
224  & *v3+4*v2+v1)+x3*y2*kyel*(4*u3+4*u2+u1)+2*x3*y3*kxel*(-4*
225  & v3-4*v2-v1)+2*x3*y3*kyel*(-4*u3-4*u2-u1)+2*y2**2*kxel*(
226  & 4*u3+4*u2+u1)+5*y2*y3*kxel*(-4*u3-4*u2-u1)+2*y3**2*kxel*
227  & (4*u3+4*u2+u1))*xmul/(54*x2*y3-54*x3*y2)
228 !
229  a34(ielem)=
230  & (2*x2**2*kyel*(-4*v3-4*v2-v1)+x2*x3*kyel*(8*v3+
231  & 11*v2-v1)+2*x2*y2*kxel*(4*v3+4*v2+v1)+2*x2*y2*kyel*(4*u3
232  & +4*u2+u1)+x2*y3*kxel*(-4*v3-4*v2-v1)+x2*y3*kyel*(-4*u3-7
233  & *u2+2*u1)+x3**2*kyel*(-8*v3-5*v2-5*v1)+x3*y2*kxel*(-4*v3
234  & -7*v2+2*v1)+x3*y2*kyel*(-4*u3-4*u2-u1)+x3*y3*kxel*(8*v3+
235  & 5*v2+5*v1)+x3*y3*kyel*(8*u3+5*u2+5*u1)+2*y2**2*kxel*(-4
236  & *u3-4*u2-u1)+y2*y3*kxel*(8*u3+11*u2-u1)+y3**2*kxel*(-8*u3
237  & -5*u2-5*u1))*xmul/(18*x2*y3-18*x3*y2)
238 !
239 ! USES HERE THE 'MAGIC SQUARE' PROPERTIES OF A DIFFUSION-LIKE MATRIX
240 ! (SUM OF EACH LINE AND EACH COLUMN IS 0)
241 !
242  a14(ielem) = - a11(ielem) - a12(ielem) - a13(ielem)
243  a22(ielem) = - a21(ielem) - a23(ielem) - a24(ielem)
244  a33(ielem) = - a31(ielem) - a32(ielem) - a34(ielem)
245  a41(ielem) = - a11(ielem) - a21(ielem) - a31(ielem)
246  a42(ielem) = - a12(ielem) - a22(ielem) - a32(ielem)
247  a43(ielem) = - a13(ielem) - a23(ielem) - a33(ielem)
248  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
249 !
250  ENDDO ! IELEM
251 !
252 !-----------------------------------------------------------------------
253 ! CASE WHERE U IS OF QUASI-BUBBLE DISCRETISATION
254 !-----------------------------------------------------------------------
255 !
256  ELSEIF(ielmf.EQ.10.AND.ielmg.EQ.10.AND.
257  & ielmu.EQ.12.AND.ielmv.EQ.12) THEN
258 !
259 ! LOOP ON THE ELEMENTS
260 !
261  DO ielem = 1 , nelem
262 !
263  x2 = xel(ielem,2)
264  x3 = xel(ielem,3)
265  y2 = yel(ielem,2)
266  y3 = yel(ielem,3)
267 !
268  u1 = u(ikle1(ielem))
269  u2 = u(ikle2(ielem))
270  u3 = u(ikle3(ielem))
271  u4 = u(ikle3(ielem))
272  v1 = v(ikle1(ielem))
273  v2 = v(ikle2(ielem))
274  v3 = v(ikle3(ielem))
275  v4 = v(ikle3(ielem))
276 !
277  kxel = f(ielem)
278  kyel = g(ielem)
279 !
280 ! COMPUTES 9 OF THE 16 TERMS
281 !
282  a11(ielem)=
283  & (x2**2*kyel*(v3+5*v4+4*v2+5*v1)+4*x2*x3*kyel*(-
284  & v3-2*v4-v2-2*v1)+x2*y2*kxel*(-v3-5*v4-4*v2-5*v1)+x2*y2
285  & *kyel*(-u3-5*u4-4*u2-5*u1)+2*x2*y3*kxel*(v3+2*v4+v2+2*
286  & v1)+2*x2*y3*kyel*(u3+2*u4+u2+2*u1)+x3**2*kyel*(4*v3+5*v4
287  & +v2+5*v1)+2*x3*y2*kxel*(v3+2*v4+v2+2*v1)+2*x3*y2*kyel*(
288  & u3+2*u4+u2+2*u1)+x3*y3*kxel*(-4*v3-5*v4-v2-5*v1)+x3*y3
289  & *kyel*(-4*u3-5*u4-u2-5*u1)+y2**2*kxel*(u3+5*u4+4*u2+5*
290  & u1)+4*y2*y3*kxel*(-u3-2*u4-u2-2*u1)+y3**2*kxel*(4*u3+5*
291  & u4+u2+5*u1))*xmul/(18*x2*y3-18*x3*y2)
292 !
293  a12(ielem)=
294  & (2*x2**2*kyel*(v4+v2+v1)+x2*x3*kyel*(v4+v2+v1)-(2
295  & *x2*y2*kxel)*(v4+v2+v1)-(2*x2*y2*kyel)*(u4+u2+u1)+x2*y3*kxel*(
296  & v4+v2+v1)-(2*x2*y3*kyel)*(u4+u2+u1)-(x3**2*kyel)*(v4+v2+v1)-
297  & (2*x3*y2*kxel)*(v4+v2+v1)+x3*y2*kyel*(u4+u2+u1)+x3*y3*kxel*(v4
298  & +v2+v1)+x3*y3*kyel*(u4+u2+u1)+2*y2**2*kxel*(u4+u2+u1)+y2*y3*
299  & kxel*(u4+u2+u1)-(y3**2*kxel)*(u4+u2+u1))
300  & *xmul/(18*x2*y3-18*x3*y2)
301 !
302  a13(ielem)=
303  & (-(x2**2*kyel)*(v3+v4+v1)+x2*x3*kyel*(v3+v4+v1)+x2*
304  & y2*kxel*(v3+v4+v1)+x2*y2*kyel*(u3+u4+u1)-(2*x2*y3*kxel)*(v3+v4
305  & +v1)+x2*y3*kyel*(u3+u4+u1)+2*x3**2*kyel*(v3+v4+v1)+x3*y2*kxel*
306  & (v3+v4+v1)-(2*x3*y2*kyel)*(u3+u4+u1)-(2*x3*y3*kxel)*(v3+v4+
307  & v1)-(2*x3*y3*kyel)*(u3+u4+u1)-(y2**2*kxel)*(u3+u4+u1)+y2*y3*
308  & kxel*(u3+u4+u1)+2*y3**2*kxel*(u3+u4+u1))
309  & *xmul/(18*x2*y3-18*x3*y2)
310 !
311  a21(ielem)=
312  & (2*x2**2*kyel*(v4+v2+v1)+x2*x3*kyel*(v4+v2+v1)-(2
313  & *x2*y2*kxel)*(v4+v2+v1)-(2*x2*y2*kyel)*(u4+u2+u1)-(2*x2*y3*
314  & kxel)*(v4+v2+v1)+x2*y3*kyel*(u4+u2+u1)-(x3**2*kyel)*(v4+v2+v1)+
315  & x3*y2*kxel*(v4+v2+v1)-(2*x3*y2*kyel)*(u4+u2+u1)+x3*y3*kxel*(v4
316  & +v2+v1)+x3*y3*kyel*(u4+u2+u1)+2*y2**2*kxel*(u4+u2+u1)+y2*y3*
317  & kxel*(u4+u2+u1)-(y3**2*kxel)*(u4+u2+u1))
318  & *xmul/(18*x2*y3-18*x3*y2)
319 !
320  a23(ielem)=
321  & (2*x2**2*kyel*(v3+v4+v2)-(5*x2*x3*kyel)*(v3+v4+v2
322  & )-(2*x2*y2*kxel)*(v3+v4+v2)-(2*x2*y2*kyel)*(u3+u4+u2)+4*x2
323  & *y3*kxel*(v3+v4+v2)+x2*y3*kyel*(u3+u4+u2)+2*x3**2*kyel*(v3+v4+
324  & v2)+x3*y2*kxel*(v3+v4+v2)+4*x3*y2*kyel*(u3+u4+u2)-(2*x3*y3*
325  & kxel)*(v3+v4+v2)-(2*x3*y3*kyel)*(u3+u4+u2)+2*y2**2*kxel*(u3+
326  & u4+u2)-(5*y2*y3*kxel)*(u3+u4+u2)+2*y3**2*kxel*(u3+u4+u2))
327  & *xmul/(18*x2*y3-18*x3*y2)
328 !
329  a24(ielem)=
330  & (x2**2*kyel*(-v3-2*v4-2*v2-v1)+x2*x3*kyel*(3*v3+
331  & 2*v4+2*v2-v1)+x2*y2*kxel*(v3+2*v4+2*v2+v1)+x2*y2*kyel*(u3+
332  & 2*u4+2*u2+u1)+x2*y3*kxel*(-2*v3-v4-v2+v1)-(x2*y3*kyel)*(u3+
333  & u4+u2)-(2*x3**2*kyel)*(v3+v4+v2)-(x3*y2*kxel)*(v3+v4+v2)+x3*
334  & y2*kyel*(-2*u3-u4-u2+u1)+2*x3*y3*kxel*(v3+v4+v2)+2*x3*y3*
335  & kyel*(u3+u4+u2)+y2**2*kxel*(-u3-2*u4-2*u2-u1)+y2*y3*kxel*(3*
336  & u3+2*u4+2*u2-u1)-(2*y3**2*kxel)*(u3+u4+u2))
337  & *xmul/(6*x2*y3-6*x3*y2)
338 !
339  a31(ielem)=
340  & (-(x2**2*kyel)*(v3+v4+v1)+x2*x3*kyel*(v3+v4+v1)+x2*
341  & y2*kxel*(v3+v4+v1)+x2*y2*kyel*(u3+u4+u1)+x2*y3*kxel*(v3+v4+v1)-
342  & (2*x2*y3*kyel)*(u3+u4+u1)+2*x3**2*kyel*(v3+v4+v1)-(2*x3*y2
343  & *kxel)*(v3+v4+v1)+x3*y2*kyel*(u3+u4+u1)-(2*x3*y3*kxel)*(v3+v4+
344  & v1)-(2*x3*y3*kyel)*(u3+u4+u1)-(y2**2*kxel)*(u3+u4+u1)+y2*y3*
345  & kxel*(u3+u4+u1)+2*y3**2*kxel*(u3+u4+u1))
346  & *xmul/(18*x2*y3-18*x3*y2)
347 !
348  a32(ielem)=
349  & (2*x2**2*kyel*(v3+v4+v2)-(5*x2*x3*kyel)*(v3+v4+v2
350  & )-(2*x2*y2*kxel)*(v3+v4+v2)-(2*x2*y2*kyel)*(u3+u4+u2)+x2*y3
351  & *kxel*(v3+v4+v2)+4*x2*y3*kyel*(u3+u4+u2)+2*x3**2*kyel*(v3+v4+
352  & v2)+4*x3*y2*kxel*(v3+v4+v2)+x3*y2*kyel*(u3+u4+u2)-(2*x3*y3*
353  & kxel)*(v3+v4+v2)-(2*x3*y3*kyel)*(u3+u4+u2)+2*y2**2*kxel*(u3+
354  & u4+u2)-(5*y2*y3*kxel)*(u3+u4+u2)+2*y3**2*kxel*(u3+u4+u2))
355  & *xmul/(18*x2*y3-18*x3*y2)
356 !
357  a34(ielem)=
358  & (-(2*x2**2*kyel)*(v3+v4+v2)+x2*x3*kyel*(2*v3+2*
359  & v4+3*v2-v1)+2*x2*y2*kxel*(v3+v4+v2)+2*x2*y2*kyel*(u3+u4+u2
360  & )-(x2*y3*kxel)*(v3+v4+v2)+x2*y3*kyel*(-u3-u4-2*u2+u1)+x3**2*
361  & kyel*(-2*v3-2*v4-v2-v1)+x3*y2*kxel*(-v3-v4-2*v2+v1)-(x3*y2
362  & *kyel)*(u3+u4+u2)+x3*y3*kxel*(2*v3+2*v4+v2+v1)+x3*y3*kyel*(2
363  & *u3+2*u4+u2+u1)-(2*y2**2*kxel)*(u3+u4+u2)+y2*y3*kxel*(2*u3
364  & +2*u4+3*u2-u1)+y3**2*kxel*(-2*u3-2*u4-u2-u1))
365  & *xmul/(6*x2*y3-6*x3*y2)
366 !
367 !
368 ! USES HERE THE 'MAGIC SQUARE' PROPERTIES OF A DIFFUSION-LIKE MATRIX
369 ! (SUM OF EACH LINE AND EACH COLUMN IS 0)
370 !
371  a14(ielem) = - a11(ielem) - a12(ielem) - a13(ielem)
372  a22(ielem) = - a21(ielem) - a23(ielem) - a24(ielem)
373  a33(ielem) = - a31(ielem) - a32(ielem) - a34(ielem)
374  a41(ielem) = - a11(ielem) - a21(ielem) - a31(ielem)
375  a42(ielem) = - a12(ielem) - a22(ielem) - a32(ielem)
376  a43(ielem) = - a13(ielem) - a23(ielem) - a33(ielem)
377  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
378 !
379  ENDDO ! IELEM
380 !
381 ! OTHER TYPES OF FUNCTIONS F AND G
382 !
383 !-----------------------------------------------------------------------
384 !
385  ELSE
386 !
387  WRITE(lu,101) ielmf,sf%NAME
388  WRITE(lu,111) ielmg,sg%NAME
389  WRITE(lu,201) ielmu,su%NAME
390  WRITE(lu,301)
391 101 FORMAT(1x,'MT03BB (BIEF) :',/,
392  & 1x,'DISCRETIZATION OF F:',1i6,
393  & 1x,'REAL NAME: ',a6)
394 111 FORMAT(1x,'DISCRETIZATION OF G:',1i6,
395  & 1x,'REAL NAME: ',a6)
396 201 FORMAT(1x,'DISCRETIZATION OF U:',1i6,
397  & 1x,'REAL NAME: ',a6)
398 301 FORMAT(1x,'CASE NOT IMPLEMENTED')
399  CALL plante(0)
400  stop
401 !
402  ENDIF
403 !
404 !-----------------------------------------------------------------------
405 !
406  RETURN
407  END
mt03bb
subroutine mt03bb(A11, A12, A13, A14, A21, A22, A23, A24, A31, A32, A33, A34, A41, A42, A43, A44, XMUL, SF, SG, SU, SV, F, G, U, V, XEL, YEL, IKLE1, IKLE2, IKLE3, NELEM, NELMAX)
Definition: mt03bb.f:12
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:36:47 for The TELEMAC-MASCARET system by doxygen 1.8.13