The TELEMAC-MASCARET system  trunk
mt04aa.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt04aa
3 ! *****************
4 !
5  &( a11 , a12 , a13 ,
6  & a22 , a23 ,
7  & a33 ,
8  & xmul,su,sv,u,v,xel,yel,surfac,ikle,nelem,nelmax)
9 !
10 !***********************************************************************
11 ! BIEF V6P1 21/08/2010
12 !***********************************************************************
13 !
14 !brief COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
15 !code
16 !+ / --> - --> --> --->
17 !+ A(I,J) = XMUL / U . GRAD(PSI1(I)) * U . GRAD(PSI2(J)) D(OMEGA)
18 !+ /OMEGA
19 !+
20 !+ U: VECTOR WITH COMPONENTS U, V
21 !
22 !warning THE JACOBIAN MUST BE POSITIVE
23 !
24 !history J-M HERVOUET (LNH)
25 !+ 12/04/93
26 !+ V5P1
27 !+
28 !
29 !history ALGIANE FROEHLY (MATMECA)
30 !+ 19/06/08
31 !+ V5P9
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 !| A22 |<--| ELEMENTS OF MATRIX
51 !| A23 |<--| ELEMENTS OF MATRIX
52 !| A33 |<--| ELEMENTS OF MATRIX
53 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
54 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
55 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
56 !| NELEM |-->| NUMBER OF ELEMENTS
57 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
58 !| SU |-->| BIEF_OBJ STRUCTURE OF U
59 !| SURFAC |-->| AREA OF TRIANGLES
60 !| SV |-->| BIEF_OBJ STRUCTURE OF V
61 !| U |-->| FUNCTION U USED IN THE FORMULA
62 !| V |-->| FUNCTION V USED IN THE FORMULA
63 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
64 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
65 !| XMUL |-->| MULTIPLICATION FACTOR
66 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
67 !
68  USE bief, ex_mt04aa => mt04aa
69 !
70  USE declarations_special
71  IMPLICIT NONE
72 !
73 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
74 !
75  INTEGER, INTENT(IN) :: NELEM,NELMAX
76  INTEGER, INTENT(IN) :: IKLE(nelmax,*)
77  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*)
78  DOUBLE PRECISION, INTENT(INOUT) :: A22(*),A23(*)
79  DOUBLE PRECISION, INTENT(INOUT) :: A33(*)
80  DOUBLE PRECISION, INTENT(IN) :: XMUL,U(*),V(*)
81  TYPE(bief_obj) , INTENT(IN) :: SU,SV
82  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
83  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
84 !
85 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
86 !
87 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
88 !
89  INTEGER IELMU,IELMV,IELEM
90 !
91  DOUBLE PRECISION SUR48,X2,X3,Y2,Y3,U1,U2,U3,U4,U5,U6
92  DOUBLE PRECISION V1,V2,V3,V4,V5,V6,AUX
93  DOUBLE PRECISION AUX1,AUX2,AUX3
94  DOUBLE PRECISION SUR144,SUR720,U123,V123,ANS1,ANS2,ANS3
95 !
96 !-----------------------------------------------------------------------
97 !
98  sur48 = xmul/48.d0
99  sur144 = xmul/144.d0
100  sur720 = xmul/720.d0
101 !
102 !-----------------------------------------------------------------------
103 !
104  ielmu = su%ELM
105  ielmv = sv%ELM
106 !
107 ! CASE WHERE U AND V ARE LINEAR
108 !
109  IF(ielmu.EQ.11.AND.ielmv.EQ.11) THEN
110 !
111 ! LOOP ON THE ELEMENTS
112 !
113  DO ielem = 1 , nelem
114 !
115  x2 = xel(ielem,2)
116  x3 = xel(ielem,3)
117  y2 = yel(ielem,2)
118  y3 = yel(ielem,3)
119 !
120  u1 = u(ikle(ielem,1))
121  u2 = u(ikle(ielem,2))
122  u3 = u(ikle(ielem,3))
123 !
124  v1 = v(ikle(ielem,1))
125  v2 = v(ikle(ielem,2))
126  v3 = v(ikle(ielem,3))
127 !
128  u123 = u1 + u2 + u3
129  v123 = v1 + v2 + v3
130 !
131  aux = sur48 / surfac(ielem)
132 !
133  aux1 = u1*u123+u2**2+u2*u3+u3**2
134  aux2 = u1*(v123+v1)+u2*(v123+v2)+u3*(v123+v3)
135  aux3 = v1*v123+v2**2+v2*v3+v3**2
136 !
137  a12(ielem) = ( 2*y3*(y2-y3) *aux1
138  & +(x3*(y3-y2)+(x3-x2)*y3)*aux2
139  & +2*x3*(x2-x3) *aux3 ) * aux
140 !
141  a13(ielem) = ( 2*y2*(y3-y2) *aux1
142  & +(2*x2*y2-x2*y3-x3*y2) *aux2
143  & +2*x2*(x3-x2) *aux3 ) * aux
144 !
145  a23(ielem) = ( -2*y2*y3 *aux1
146  & +(x2*y3+x3*y2) *aux2
147  & -2*x2*x3 *aux3 ) * aux
148 !
149 ! USES HERE THE 'MAGIC SQUARE' PROPERTIES
150 ! (SUM OF EACH LINE = 0 FOR EXAMPLE)
151 ! AND THE SYMMETRICAL PROPERTIES OF THE MATRIX
152 !
153  a11(ielem) = - a12(ielem) - a13(ielem)
154  a22(ielem) = - a23(ielem) - a12(ielem)
155  a33(ielem) = - a13(ielem) - a23(ielem)
156 !
157  ENDDO ! IELEM
158 !
159 !-----------------------------------------------------------------------
160 ! CASE WHERE U AND V ARE QUASI-BUBBLE
161 !
162  ELSEIF(ielmu.EQ.12.AND.ielmv.EQ.12) THEN
163 !
164 ! LOOP ON THE ELEMENTS
165 !
166  DO ielem = 1 , nelem
167 !
168  x2 = xel(ielem,2)
169  x3 = xel(ielem,3)
170  y2 = yel(ielem,2)
171  y3 = yel(ielem,3)
172 !
173  u1 = u(ikle(ielem,1))
174  u2 = u(ikle(ielem,2))
175  u3 = u(ikle(ielem,3))
176  u4 = u(ikle(ielem,4))
177 !
178  v1 = v(ikle(ielem,1))
179  v2 = v(ikle(ielem,2))
180  v3 = v(ikle(ielem,3))
181  v4 = v(ikle(ielem,4))
182 !
183  aux = sur144 / surfac(ielem)
184 !
185  a12(ielem) = (x2*((2*((v4+v2)*v3+v3**2+v4**2+v4*v2+v2**2)*
186  & x3-(2*u3+u4+u2)*v3*y3-(u3+2*u4+u2)*v4*y3-(u3+u4+2*u2)*
187  & v2*y3)+(2*((v4+v1)*v3+v3**2+v4**2+v4*v1+v1**2)*x3-(2*u3
188  & +u4+u1)*v3*y3-(u3+2*u4+u1)*v4*y3-(u3+u4+2*u1)*v1*y3)+(
189  & 2*((v2+v1)*v4+v4**2+v2**2+v2*v1+v1**2)*x3-(2*u4+u2+u1)*
190  & v4*y3-(u4+2*u2+u1)*v2*y3-(u4+u2+2*u1)*v1*y3))-(2*x3**2
191  & )*(((v4+v2)*v3+v3**2+v4**2+v4*v2+v2**2)+((v4+v1)*v3+v3**2
192  & +v4**2+v4*v1+v1**2)+((v2+v1)*v4+v4**2+v2**2+v2*v1+v1**2))
193  & +x3*(2*y3-y2)*(((2*u3+u4+u2)*v3+(u3+2*u4+u2)*v4+(u3+u4
194  & +2*u2)*v2)+((2*u3+u4+u1)*v3+(u3+2*u4+u1)*v4+(u3+u4+2*
195  & u1)*v1)+((2*u4+u2+u1)*v4+(u4+2*u2+u1)*v2+(u4+u2+2*u1)*
196  & v1))+2*y3*(y3-y2)*(-(u4+u2)*u3-(u4+u1)*u3-(u2+u1)*u4-2*
197  & u3**2-3*u4**2-u4*u2-u4*u1-2*u2**2-u2*u1-2*u1**2))*aux
198 !
199  a13(ielem) = (-(2*x2**2)*(((v4+v2)*v3+v3**2+v4**2+v4*v2+v2
200  & **2)+((v4+v1)*v3+v3**2+v4**2+v4*v1+v1**2)+((v2+v1)*v4+v4
201  & **2+v2**2+v2*v1+v1**2))+x2*((2*((v4+v2)*v3+v3**2+v4**2+
202  & v4*v2+v2**2)*x3-(2*u3+u4+u2)*(y3-2*y2)*v3-(u3+2*u4+u2)
203  & *(y3-2*y2)*v4-(u3+u4+2*u2)*(y3-2*y2)*v2)+(2*((v4+v1)*
204  & v3+v3**2+v4**2+v4*v1+v1**2)*x3-(2*u3+u4+u1)*(y3-2*y2)*
205  & v3-(u3+2*u4+u1)*(y3-2*y2)*v4-(u3+u4+2*u1)*(y3-2*y2)*
206  & v1)+(2*((v2+v1)*v4+v4**2+v2**2+v2*v1+v1**2)*x3-(2*u4+u2
207  & +u1)*(y3-2*y2)*v4-(u4+2*u2+u1)*(y3-2*y2)*v2-(u4+u2+2*
208  & u1)*(y3-2*y2)*v1))-(x3*y2)*(((2*u3+u4+u2)*v3+(u3+2*u4+
209  & u2)*v4+(u3+u4+2*u2)*v2)+((2*u3+u4+u1)*v3+(u3+2*u4+u1)*
210  & v4+(u3+u4+2*u1)*v1)+((2*u4+u2+u1)*v4+(u4+2*u2+u1)*v2+(
211  & u4+u2+2*u1)*v1))+2*y2*(y3-y2)*((u4+u2)*u3+(u4+u1)*u3+(
212  & u2+u1)*u4+2*u3**2+3*u4**2+u4*u2+u4*u1+2*u2**2+u2*u1+2
213  & *u1**2))*aux
214 !
215  a23(ielem) = (2*x2*x3*(-2*v3**2-2*v3*v4-v3*v2-v3*v1-3*v4
216  & **2-2*v4*v2-2*v4*v1-2*v2**2-v2*v1-2*v1**2)+x2*y3*(4*
217  & v3*u3+2*v3*u4+v3*u2+v3*u1+2*v4*u3+6*v4*u4+2*v4*u2+2*
218  & v4*u1+v2*u3+2*v2*u4+4*v2*u2+v2*u1+v1*u3+2*v1*u4+v1*u2+
219  & 4*v1*u1)+x3*y2*(4*v3*u3+2*v3*u4+v3*u2+v3*u1+2*v4*u3+6
220  & *v4*u4+2*v4*u2+2*v4*u1+v2*u3+2*v2*u4+4*v2*u2+v2*u1+v1
221  & *u3+2*v1*u4+v1*u2+4*v1*u1)+2*y2*y3*(-2*u3**2-2*u3*u4
222  & -u3*u2-u3*u1-3*u4**2-2*u4*u2-2*u4*u1-2*u2**2-u2*u1-2
223  & *u1**2))*aux
224 !
225 ! USES HERE THE 'MAGIC SQUARE' PROPERTIES
226 ! (SUM OF EACH LINE = 0 FOR EXAMPLE)
227 ! AND THE SYMMETRICAL PROPERTIES OF THE MATRIX
228 !
229  a11(ielem) = - a12(ielem) - a13(ielem)
230  a22(ielem) = - a23(ielem) - a12(ielem)
231  a33(ielem) = - a13(ielem) - a23(ielem)
232 !
233  ENDDO ! IELEM
234 !
235 !
236 !-----------------------------------------------------------------------
237 ! CASE WHERE U AND V ARE P2
238 !
239  ELSEIF(ielmu.EQ.13.AND.ielmv.EQ.13) THEN
240 !
241 ! LOOP ON THE ELEMENTS
242 !
243  DO ielem = 1 , nelem
244 !
245  x2 = xel(ielem,2)
246  x3 = xel(ielem,3)
247  y2 = yel(ielem,2)
248  y3 = yel(ielem,3)
249 !
250  u1 = u(ikle(ielem,1))
251  u2 = u(ikle(ielem,2))
252  u3 = u(ikle(ielem,3))
253  u4 = u(ikle(ielem,4))
254  u5 = u(ikle(ielem,5))
255  u6 = u(ikle(ielem,6))
256 !
257  v1 = v(ikle(ielem,1))
258  v2 = v(ikle(ielem,2))
259  v3 = v(ikle(ielem,3))
260  v4 = v(ikle(ielem,4))
261  v5 = v(ikle(ielem,5))
262  v6 = v(ikle(ielem,6))
263 !
264  aux = sur720 / surfac(ielem)
265 !
266  ans1 = 3.d0*u2**2*y2**2-4.d0*u4*y3**2*u3+v2*x2*u3*y2-
267  & 4.d0*v1*x3**2*v5+u2*y3*v3*x3-u2*y3**2*u3-4.d0*v4*x3**2*v3-
268  & v2*x2**2*v3-4.d0*v2*x2**2*v6-v2*x2**2*v1-v2*x3**2*v1-
269  & 4.d0*u4*y2**2*u3-v1*x2**2*v3-6.d0*u3**2*y3*y2-
270  & 32.d0*u5**2*y2*y3+16.d0*u4*y3**2*u6+16.d0*v6*x3**2*v5+
271  & 16.d0*u6*y3**2*u5+16.d0*v6**2*x3**2+16.d0*u4*y2**2*u5-
272  & 4.d0*u1*y3**2*u5-32.d0*u6**2*y2*y3+3.d0*v1**2*x3**2+
273  & 16.d0*v4*x2**2*v5-4.d0*v1*x2**2*v5+16.d0*v5**2*x3**2-
274  & 6.d0*u1**2*y3*y2-u1*y3**2*u3-32.d0*v6**2*x3*x2-
275  & 32.d0*u4**2*y2*y3+16.d0*u4*y2**2*u6+16.d0*v4**2*x2**2-
276  & 4.d0*u2*y3**2*u6+16.d0*u5**2*y3**2+3.d0*v2**2*x3**2-
277  & 6.d0*v1**2*x2*x3+3.d0*u3**2*y3**2+8.d0*u2*y2*u6*y3+
278  & 2.d0*v2*x2*v3*x3+2.d0*u2*y2*u3*y3-u2*y2*v1*x3-
279  & 6.d0*v3**2*x2*x3+16.d0*v4*x3**2*v6-4.d0*u2*y2*v6*x3+
280  & 3.d0*v3**2*x3**2+16.d0*u4**2*y3**2-6.d0*u2**2*y2*y3-
281  & 6.d0*v2**2*x2*x3+16.d0*v4*x2**2*v6-4.d0*v4*x2**2*v3+
282  & 6.d0*u2*y2*v2*x3+6.d0*v2*x2*u2*y3+32.d0*u5*y2*v5*x3+
283  & 8.d0*u4*y3*u3*y2-6.d0*u1*y3*v1*x3+4.d0*v4*x3*u3*y3-
284  & 16.d0*v4*x3*u6*y3+u1*y3*v3*x3-16.d0*u4*y3*v5*x3
285  ans2 = -32.d0*v4*x3*u4*y3-16.d0*v4*x3*u5*y3-6.d0*v3*x3*u3*y3-
286  & 32.d0*u5*y3*v5*x3+4.d0*u4*y3*v3*x3-16.d0*u4*y3*v6*x3-
287  & 16.d0*u6*y3*v5*x3+4.d0*u1*y3*v5*x3+4.d0*v1*x3*u5*y3-
288  & 32.d0*v6*x3*u6*y3-16.d0*v6*x3*u5*y3+v1*x3*u3*y3+
289  & 4.d0*v2*x3*u6*y3+4.d0*u2*y3*v6*x3+u2*y3*v1*x3+v2*x3*u1*y3+
290  & v2*x3*u3*y3+3.d0*u1**2*y3**2+3.d0*u3**2*y2**2+
291  & 3.d0*v3**2*x2**2-6.d0*u2*y3*v2*x3-4.d0*u4*y3*v3*x2-
292  & 32.d0*u4*y3*u6*y2+32.d0*v6*x3*u6*y2-32.d0*v5**2*x2*x3+
293  & 3.d0*u2**2*y3**2+16.d0*v4**2*x3**2+16.d0*u6*y3*v5*x2+
294  & 8.d0*u1*y3*u5*y2-4.d0*u1*y3*v5*x2-32.d0*u6*y2*u5*y3-
295  & 4.d0*u1*y2*v5*x3+16.d0*u6*y2*v5*x3-4.d0*v1*x3*u5*y2+
296  & 32.d0*v6*x2*u6*y3+16.d0*v6*x2*u5*y3-32.d0*v6*x3*v5*x2-
297  & 4.d0*v1*x2*u5*y3+8.d0*v1*x2*v5*x3+16.d0*v6*x3*u5*y2-
298  & v1*x3*u3*y2+6.d0*u1*y2*v1*x3-v2*x3*u1*y2-u1*y2*v3*x3-
299  & 4.d0*v2*x3*u6*y2-u2*y3*v1*x2-4.d0*u2*y3*v6*x2-v2*x3*u3*y2-
300  & 4.d0*u4*y2*v3*x3+16.d0*u4*y2*v6*x3-u2*y3*v3*x2-
301  & v1*x2*u3*y3+8.d0*v2*x2*v6*x3+16.d0*u6**2*y3**2-
302  & v2*x2*u3*y3-u2*y2*v3*x3+2.d0*u2*y2*u1*y3-v2*x2*u1*y3+
303  & 2.d0*v2*x2*v1*x3-4.d0*v2*x2*u6*y3-32.d0*v4*x2*v5*x3
304  ans3 = 16.d0*v4*x2*u5*y3+32.d0*v4*x2*u4*y3+2.d0*u1*y3*u3*y2-
305  & u1*y3*v3*x2+16.d0*u4*y3*v6*x2+16.d0*v4*x2*u6*y3+
306  & 16.d0*v4*x3*u6*y2+8.d0*v4*x2*v3*x3-32.d0*v4*x2*v6*x3-
307  & 4.d0*v4*x2*u3*y3+2.d0*v1*x2*v3*x3-4.d0*v4*x3*u3*y2+
308  & 6.d0*u1*y3*v1*x2+6.d0*v3*x2*u3*y3+16.d0*u4*y3*v5*x2+
309  & 16.d0*v4*x3*u5*y2+16.d0*u6*y2**2*u5-32.d0*u4*y2*u5*y3+
310  & 32.d0*u4*y2*v4*x3+16.d0*u4*y2*v5*x3+16.d0*u4*y3**2*u5-
311  & 6.d0*u2*y2*v2*x2+u2*y2*v1*x2+4.d0*u2*y2*v6*x2+
312  & 6.d0*v3*x3*u3*y2+32.d0*u5*y3*v5*x2-v1*x3**2*v3
313  & -32.d0*u5*y2*v5*x2+16.d0*v6*x2**2*v5-u1*y2**2*u3-
314  & u2*y2**2*u1-4.d0*u1*y2**2*u5+3.d0*v2**2*x2**2-
315  & v2*x3**2*v3+16.d0*u6**2*y2**2+4.d0*u1*y2*v5*x2-
316  & 16.d0*u6*y2*v5*x2-4.d0*u2*y2**2*u6+3.d0*v1**2*x2**2+
317  & 16.d0*v6**2*x2**2-32.d0*v4**2*x2*x3+16.d0*v4*x3**2*v5+
318  & v1*x2*u3*y2-6.d0*v1*x2*u1*y2+4.d0*v4*x2*u3*y2-
319  & 16.d0*v4*x2*u6*y2+u2*y2*v3*x2+4.d0*v2*x2*u6*y2+
320  & v2*x2*u1*y2-16.d0*v4*x2*u5*y2+3.d0*u1**2*y2**2+
321  & u1*y2*v3*x2+4.d0*u4*y2*v3*x2-16.d0*u4*y2*v6*x2-
322  & 32.d0*v6*x2*u6*y2-16.d0*v6*x2*u5*y2+4.d0*v1*x2*u5*y2
323  a11(ielem) = (16.d0*u4**2*y2**2-u2*y3**2*u1-4*v2*x3**2*v6-
324  & u2*y2**2*u3+16.d0*v5**2*x2**2-32.d0*u4*y2*v4*x2-
325  & 6.d0*v3*x2*u3*y2-16.d0*u4*y2*v5*x2+
326  & 16.d0*u5**2*y2**2 + ans1 + ans2 + ans3)*2.d0*aux
327 !
328  a22(ielem) = (-4.d0*u4*y3**2*u3-4.d0*v1*x3**2*v5+u2*y3*v3*x3-
329  & u2*y3**2*u3-4.d0*v4*x3**2*v3-v2*x3**2*v1+
330  & 16.d0*u4*y3**2*u6+16.d0*v6*x3**2*v5+16.d0*u6*y3**2*u5+
331  & 16.d0*v6**2*x3**2-4.d0*u1*y3**2*u5+3*v1**2*x3**2+
332  & 16.d0*v5**2*x3**2-u1*y3**2*u3-4.d0*u2*y3**2*u6+
333  & 16.d0*u5**2*y3**2+3.d0*v2**2*x3**2+3.d0*u3**2*y3**2+
334  & 16.d0*v4*x3**2*v6+3.d0*v3**2*x3**2+16.d0*u4**2*y3**2-
335  & 6.d0*u1*y3*v1*x3+4.d0*v4*x3*u3*y3-16.d0*v4*x3*u6*y3+
336  & u1*y3*v3*x3-16.d0*u4*y3*v5*x3-32.d0*v4*x3*u4*y3-
337  & 16.d0*v4*x3*u5*y3-6.d0*v3*x3*u3*y3-32.d0*u5*y3*v5*x3+
338  & 4.d0*u4*y3*v3*x3-16.d0*u4*y3*v6*x3-16.d0*u6*y3*v5*x3+
339  & 4.d0*u1*y3*v5*x3+4.d0*v1*x3*u5*y3-32.d0*v6*x3*u6*y3-
340  & 16.d0*v6*x3*u5*y3+v1*x3*u3*y3+4.d0*v2*x3*u6*y3+
341  & 4.d0*u2*y3*v6*x3+u2*y3*v1*x3+v2*x3*u1*y3+v2*x3*u3*y3+
342  & 3.d0*u1**2*y3**2-6.d0*u2*y3*v2*x3+3.d0*u2**2*y3**2+
343  & 16.d0*v4**2*x3**2+16.d0*u6**2*y3**2+16.d0*u4*y3**2*u5-
344  & v1*x3**2*v3-v2*x3**2*v3+16.d0*v4*x3**2*v5-u2*y3**2*u1-
345  & 4.d0*v2*x3**2*v6)*2.d0*aux
346 !
347  ans1 = 6.d0*u3**2*y3*y2+32.d0*u5**2*y2*y3+32.d0*u6**2*y2*y3+
348  & 6.d0*u1**2*y3*y2+32.d0*v6**2*x3*x2+32.d0*u4**2*y2*y3+
349  & 6.d0*v1**2*x2*x3-8.d0*u2*y2*u6*y3-2.d0*v2*x2*v3*x3-
350  & 2.d0*u2*y2*u3*y3+u2*y2*v1*x3+6.d0*v3**2*x2*x3+
351  & 4.d0*u2*y2*v6*x3+6.d0*u2**2*y2*y3+6.d0*v2**2*x2*x3-
352  & 6.d0*u2*y2*v2*x3-6.d0*v2*x2*u2*y3-32.d0*u5*y2*v5*x3-
353  & 8.d0*u4*y3*u3*y2+4.d0*u4*y3*v3*x2+32.d0*u4*y3*u6*y2-
354  & 32.d0*v6*x3*u6*y2+32.d0*v5**2*x2*x3-16.d0*u6*y3*v5*x2-
355  & 8.d0*u1*y3*u5*y2+4.d0*u1*y3*v5*x2+32.d0*u6*y2*u5*y3+
356  & 4.d0*u1*y2*v5*x3-16.d0*u6*y2*v5*x3+4.d0*v1*x3*u5*y2-
357  & 32.d0*v6*x2*u6*y3-16.d0*v6*x2*u5*y3+32.d0*v6*x3*v5*x2+
358  & 4.d0*v1*x2*u5*y3-8.d0*v1*x2*v5*x3-16.d0*v6*x3*u5*y2+
359  & v1*x3*u3*y2-6.d0*u1*y2*v1*x3+v2*x3*u1*y2+u1*y2*v3*x3+
360  & 4.d0*v2*x3*u6*y2+u2*y3*v1*x2+4.d0*u2*y3*v6*x2+v2*x3*u3*y2+
361  & 4.d0*u4*y2*v3*x3-16.d0*u4*y2*v6*x3+u2*y3*v3*x2+v1*x2*u3*y3-
362  & 8.d0*v2*x2*v6*x3+v2*x2*u3*y3+u2*y2*v3*x3-2.d0*u2*y2*u1*y3+
363  & v2*x2*u1*y3-2.d0*v2*x2*v1*x3+4.d0*v2*x2*u6*y3+
364  & 32.d0*v4*x2*v5*x3-16.d0*v4*x2*u5*y3-32.d0*v4*x2*u4*y3-
365  & 2.d0*u1*y3*u3*y2+u1*y3*v3*x2-16.d0*u4*y3*v6*x2
366  a23(ielem) = -(-16.d0*v4*x2*u6*y3-16.d0*v4*x3*u6*y2-
367  & 8.d0*v4*x2*v3*x3+32.d0*v4*x2*v6*x3+
368  & 4.d0*v4*x2*u3*y3-2.d0*v1*x2*v3*x3+
369  & 4.d0*v4*x3*u3*y2-6.d0*u1*y3*v1*x2-
370  & 6.d0*v3*x2*u3*y3-16.d0*u4*y3*v5*x2-
371  & 16.d0*v4*x3*u5*y2+32.d0*u4*y2*u5*y3-
372  & 32.d0*u4*y2*v4*x3-16.d0*u4*y2*v5*x3-
373  & 6.d0*v3*x3*u3*y2-32.d0*u5*y3*v5*x2+
374  & 32.d0*v4**2*x2*x3 + ans1)*aux
375 !
376 ! USES HERE THE 'MAGIC SQUARE' PROPERTIES
377 ! (SUM OF EACH LINE = 0 FOR EXAMPLE)
378 ! AND THE SYMMETRICAL PROPERTIES OF THE MATRIX
379 !
380  a12(ielem) = - a22(ielem) - a23(ielem)
381  a13(ielem) = - a11(ielem) - a12(ielem)
382  a33(ielem) = - a13(ielem) - a23(ielem)
383 !
384  ENDDO ! IELEM
385 !
386 !
387 ! OTHER TYPES OF DISCRETISATION FOR U
388 !
389 !-----------------------------------------------------------------------
390 !
391  ELSE
392 !
393  IF(ielmu.EQ.ielmv) THEN
394  WRITE(lu,101) ielmu
395 101 FORMAT(1x,'MT04AA (BIEF) :',/,
396  & 1x,'DISCRETIZATION OF U AND V : ',1i6,' NOT AVAILABLE')
397  ELSE
398  WRITE(lu,201) ielmu,ielmv
399 201 FORMAT(1x,'MT04AA (BIEF) :',/,
400  & 1x,'U AND V OF A DIFFERENT DISCRETISATION:',1i6,3x,1i6)
401  ENDIF
402 !
403  CALL plante(1)
404  stop
405 !
406  ENDIF
407 !
408 !-----------------------------------------------------------------------
409 !
410  RETURN
411  END
declarations_special
Definition: declarations_special.F:3
mt04aa
subroutine mt04aa(A11, A12, A13, A22, A23, A33, XMUL, SU, SV, U, V, XEL, YEL, SURFAC, IKLE, NELEM, NELMAX)
Definition: mt04aa.f:10
declarations_special::lu
integer lu
Definition: declarations_special.F:45
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:36:53 for The TELEMAC-MASCARET system by doxygen 1.8.13