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