The TELEMAC-MASCARET system  trunk
mt12ba.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt12ba
3 ! *****************
4 !
5  &( a11 , a12 , a13 ,
6  & a21 , a22 , a23 ,
7  & a31 , a32 , a33 ,
8  & a41 , a42 , a43 ,
9  & xmul,sf,su,sv,f,u,v,
10  & xel,yel,surfac,
11  & ikle1,ikle2,ikle3,ikle4,
12  & nelem,nelmax,icoord)
13 !
14 !***********************************************************************
15 ! BIEF V6P1 21/08/2010
16 !***********************************************************************
17 !
18 !brief COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
19 !code
20 !+ EXAMPLE WITH ICOORD=1
21 !+
22 !+ / D -> -->
23 !+ A(I,J)= XMUL / PSI2(J) * --(F) * U.GRAD(PSI1(I)) D(OMEGA)
24 !+ /OMEGA DX
25 !+
26 !+
27 !+ PSI1 : BASES OF TYPE QUASI-BUBBLE
28 !+ PSI2 : BASES OF TYPE LINEAR
29 !+ F : FUNCTION OF TYPE QUASI-BUBBLE TRIANGLE
30 !+ U : VECTOR OF TYPE LINEAR OR P0
31 !+
32 !+ IT WOULD BE A DERIVATIVE WRT Y WITH ICOORD=2
33 !
34 !warning THE JACOBIAN MUST BE POSITIVE
35 !
36 !history J-M HERVOUET (LNH) ; C MOULIN (LNH)
37 !+ 09/12/94
38 !+ V5P1
39 !+
40 !
41 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
42 !+ 13/07/2010
43 !+ V6P0
44 !+ Translation of French comments within the FORTRAN sources into
45 !+ English comments
46 !
47 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
48 !+ 21/08/2010
49 !+ V6P0
50 !+ Creation of DOXYGEN tags for automated documentation and
51 !+ cross-referencing of the FORTRAN sources
52 !
53 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
54 !| A11 |<--| ELEMENTS OF MATRIX
55 !| A12 |<--| ELEMENTS OF MATRIX
56 !| A13 |<--| ELEMENTS OF MATRIX
57 !| A21 |<--| ELEMENTS OF MATRIX
58 !| A22 |<--| ELEMENTS OF MATRIX
59 !| A23 |<--| ELEMENTS OF MATRIX
60 !| A31 |<--| ELEMENTS OF MATRIX
61 !| A32 |<--| ELEMENTS OF MATRIX
62 !| A33 |<--| ELEMENTS OF MATRIX
63 !| A41 |<--| ELEMENTS OF MATRIX
64 !| A42 |<--| ELEMENTS OF MATRIX
65 !| A43 |<--| 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 !| IKLE4 |-->| FOURTH POINTS OF TRIANGLES (QUADRATIC)
72 !| IKLE5 |-->| FIFTH POINTS OF TRIANGLES (QUADRATIC)
73 !| IKLE6 |-->| SIXTH POINTS OF TRIANGLES (QUADRATIC)
74 !| NELEM |-->| NUMBER OF ELEMENTS
75 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
76 !| SF |-->| STRUCTURE OF FUNCTIONS F
77 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
78 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
79 !| XMUL |-->| MULTIPLICATION FACTOR
80 !| F |-->| FUNCTION USED IN THE FORMULA
81 !| ICOORD |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
82 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
83 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
84 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
85 !| IKLE4 |-->| QUASI-BUBBLE POINT
86 !| NELEM |-->| NUMBER OF ELEMENTS
87 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
88 !| SF |-->| STRUCTURE OF FUNCTIONS F
89 !| SU |-->| BIEF_OBJ STRUCTURE OF U
90 !| SURFAC |-->| AREA OF TRIANGLES
91 !| SV |-->| BIEF_OBJ STRUCTURE OF V
92 !| U |-->| FUNCTION U USED IN THE FORMULA
93 !| V |-->| FUNCTION V USED IN THE FORMULA
94 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
95 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
96 !| XMUL |-->| MULTIPLICATION FACTOR
97 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
98 !
99  USE bief, ex_mt12ba => mt12ba
100 !
101  USE declarations_special
102  IMPLICIT NONE
103 !
104 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
105 !
106  INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
107  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax)
108  INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
109 !
110  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*)
111  DOUBLE PRECISION, INTENT(INOUT) :: A21(*),A22(*),A23(*)
112  DOUBLE PRECISION, INTENT(INOUT) :: A31(*),A32(*),A33(*)
113  DOUBLE PRECISION, INTENT(INOUT) :: A41(*),A42(*),A43(*)
114 !
115  DOUBLE PRECISION, INTENT(IN) :: XMUL
116  DOUBLE PRECISION, INTENT(IN) :: F(*),U(*),V(*)
117 !
118 ! STRUCTURES OF F, U, V
119  TYPE(bief_obj), INTENT(IN) :: SF,SU,SV
120 !
121 !
122  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
123  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
124 !
125 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
126 !
127  INTEGER IELEM,IELMF,IELMU,IELMV
128  DOUBLE PRECISION X2,X3,Y2,Y3,F1,F2,F3,F4
129  DOUBLE PRECISION U1,U2,U3,V1,V2,V3,UX,UY,AUX108,XSU108,AUX036
130  DOUBLE PRECISION AX1296,XS1296,XSUR36,AUX432,XSU432
131 !
132 !-----------------------------------------------------------------------
133 !
134  ielmf=sf%ELM
135  ielmu=su%ELM
136  ielmv=sv%ELM
137 !
138  xsur36 = xmul / 36.d0
139  xsu108 = xmul / 108.d0
140  xsu432 = xmul / 432.d0
141  xs1296 = xmul /1296.d0
142 !
143 !-----------------------------------------------------------------------
144 ! CASE WHERE F IS OF TYPE P1
145 !-----------------------------------------------------------------------
146 !
147  IF(ielmf.EQ.12) THEN
148 !
149  IF(ielmu.EQ.10.AND.ielmv.EQ.10) THEN
150 !
151 !================================
152 ! DERIVATIVE WRT X =
153 !================================
154 !
155  IF(icoord.EQ.1) THEN
156 !
157 ! LOOP ON THE ELEMENTS
158 !
159  DO ielem = 1 , nelem
160 !
161 ! INITIALISES THE GEOMETRICAL VARIABLES
162 !
163  x2 = xel(ielem,2)
164  x3 = xel(ielem,3)
165  y2 = yel(ielem,2)
166  y3 = yel(ielem,3)
167 !
168  f1 = f(ikle1(ielem))
169  f2 = f(ikle2(ielem)) - f1
170  f3 = f(ikle3(ielem)) - f1
171  f4 = f(ikle4(ielem)) - f1
172 !
173  ux = u(ielem)
174  uy = v(ielem)
175 !
176  aux108 = xsu108 / surfac(ielem)
177  aux036 = xsur36 / surfac(ielem)
178 !
179 ! EXTRADIAGONAL TERMS
180 !
181  a12(ielem)=((((f3-3*f4)*y3+y2*f3)*x2*uy-2*((f3-3*f4)*y3
182  & +y2*f3)*x3*uy+8*((3*f4-f2)*y2-y3*f2)*x2*uy-4*((3*f4-
183  & f2)*y2-y3*f2)*x3*uy+(y3*f3-3*y3*f4+y2*f3)*(2*y3-y2)*ux-
184  & 4*(y3*f2-3*y2*f4+y2*f2)*(y3-2*y2)*ux))*aux108
185  a13(ielem)=((4*((f3-3*f4)*y3+y2*f3)*x2*uy-8*((f3-3*f4)
186  & *y3+y2*f3)*x3*uy+2*((3*f4-f2)*y2-y3*f2)*x2*uy-((3*f4-
187  & f2)*y2-y3*f2)*x3*uy+4*(y3*f3-3*y3*f4+y2*f3)*(2*y3-y2)*
188  & ux-(y3*f2-3*y2*f4+y2*f2)*(y3-2*y2)*ux))*aux108
189  a21(ielem)=(-(((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*x2
190  & *uy-2*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*x3*uy-4
191  & *((3*f4-f2)*y2-y3*f2)*x2*uy-4*((3*f4-f2)*y2-y3*f2)*x3*
192  & uy-(y3*f3-3*y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-y2*f2)*(2*
193  & y3-y2)*ux-4*(y3*f2-3*y2*f4+y2*f2)*(y3+y2)*ux))*aux108
194  a23(ielem)=(-(4*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)
195  & *x2*uy-8*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*x3*uy
196  & -((3*f4-f2)*y2-y3*f2)*x2*uy-((3*f4-f2)*y2-y3*f2)*x3*uy-
197  & 4*(y3*f3-3*y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-y2*f2)*(2*
198  & y3-y2)*ux-(y3*f2-3*y2*f4+y2*f2)*(y3+y2)*ux))*aux108
199  a31(ielem)=(-(2*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)
200  & *x2*uy-((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*x3*uy+4
201  & *((f3-3*f4)*y3+y2*f3)*x2*uy+4*((f3-3*f4)*y3+y2*f3)*x3*
202  & uy-(y3*f3-3*y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-y2*f2)*(y3-
203  & 2*y2)*ux-4*(y3*f3-3*y3*f4+y2*f3)*(y3+y2)*ux))*aux108
204  a32(ielem)=(-(8*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)
205  & *x2*uy-4*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*x3*uy
206  & +((f3-3*f4)*y3+y2*f3)*x2*uy+((f3-3*f4)*y3+y2*f3)*x3*uy-
207  & 4*(y3*f3-3*y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-y2*f2)*(y3-
208  & 2*y2)*ux-(y3*f3-3*y3*f4+y2*f3)*(y3+y2)*ux))*aux108
209  a41(ielem)=(-(x2*uy*y3*f3-3*x2*uy*y3*f4-2*x2*uy*y3*f2-2
210  & *x2*uy*y2*f3+15*x2*uy*y2*f4-5*x2*uy*y2*f2-5*x3*uy*y3*
211  & f3+15*x3*uy*y3*f4-2*x3*uy*y3*f2-2*x3*uy*y2*f3-3*x3*uy
212  & *y2*f4+x3*uy*y2*f2+5*ux*y3**2*f3-15*ux*y3**2*f4+2*ux*
213  & y3**2*f2+ux*y3*y2*f3+6*ux*y3*y2*f4+ux*y3*y2*f2+2*ux*y2
214  & **2*f3-15*ux*y2**2*f4+5*ux*y2**2*f2))*aux036
215  a42(ielem)=((4*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*
216  & x2*uy-4*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*x3*uy+
217  & ((f3-3*f4)*y3+y2*f3)*x3*uy-((f3-3*f4)*y3+y2*f3)*ux*y3-
218  & 4*((3*f4-f2)*y2-y3*f2)*x2*uy+4*((3*f4-f2)*y2-y3*f2)*ux
219  & *y2-4*(y3*f3-3*y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-y2*f2)*
220  & (y3-y2)*ux))*aux036
221  a43(ielem)=((4*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*
222  & x2*uy-4*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*x3*uy+
223  & 4*((f3-3*f4)*y3+y2*f3)*x3*uy-4*((f3-3*f4)*y3+y2*f3)*ux
224  & *y3-((3*f4-f2)*y2-y3*f2)*x2*uy+((3*f4-f2)*y2-y3*f2)*ux*
225  & y2-4*(y3*f3-3*y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-y2*f2)*(
226  & y3-y2)*ux))*aux036
227 !
228 ! DIAGONAL TERMS
229 ! THE SUM OF EACH COLUMN IS 0
230 !
231  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
232  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
233  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
234 !
235  ENDDO ! IELEM
236 !
237  ELSEIF(icoord.EQ.2) THEN
238 !
239 !================================
240 ! DERIVATIVE WRT Y =
241 !================================
242 !
243  DO ielem = 1 , nelem
244 !
245 ! INITIALISES THE GEOMETRICAL VARIABLES
246 !
247  x2 = xel(ielem,2)
248  x3 = xel(ielem,3)
249  y2 = yel(ielem,2)
250  y3 = yel(ielem,3)
251 !
252  f1 = f(ikle1(ielem))
253  f2 = f(ikle2(ielem)) - f1
254  f3 = f(ikle3(ielem)) - f1
255  f4 = f(ikle4(ielem)) - f1
256 !
257  ux = u(ielem)
258  uy = v(ielem)
259 !
260  aux108 = xsu108 / surfac(ielem)
261  aux036 = xsur36 / surfac(ielem)
262 !
263 ! EXTRADIAGONAL TERMS
264 !
265  a12(ielem)=(-(x2**2*uy*f3+24*x2**2*uy*f4-8*x2**2*uy*f2-
266  & x2*x3*uy*f3-15*x2*x3*uy*f4-4*x2*x3*uy*f2+2*x2*ux*y3*f3
267  & +12*x2*ux*y3*f4-4*x2*ux*y3*f2-x2*ux*y2*f3-24*x2*ux*y2*
268  & f4+8*x2*ux*y2*f2-2*x3**2*uy*f3+6*x3**2*uy*f4+4*x3**2*
269  & uy*f2+2*x3*ux*y3*f3-6*x3*ux*y3*f4-4*x3*ux*y3*f2-x3*ux*
270  & y2*f3+3*x3*ux*y2*f4+8*x3*ux*y2*f2))*aux108
271  a13(ielem)=(-(4*x2**2*uy*f3+6*x2**2*uy*f4-2*x2**2*uy*f2
272  & -4*x2*x3*uy*f3-15*x2*x3*uy*f4-x2*x3*uy*f2+8*x2*ux*y3*
273  & f3+3*x2*ux*y3*f4-x2*ux*y3*f2-4*x2*ux*y2*f3-6*x2*ux*y2*
274  & f4+2*x2*ux*y2*f2-8*x3**2*uy*f3+24*x3**2*uy*f4+x3**2*uy
275  & *f2+8*x3*ux*y3*f3-24*x3*ux*y3*f4-x3*ux*y3*f2-4*x3*ux*
276  & y2*f3+12*x3*ux*y2*f4+2*x3*ux*y2*f2))*aux108
277  a21(ielem)=((2*x2**2*uy*f3-15*x2**2*uy*f4+5*x2**2*uy*f2
278  & -5*x2*x3*uy*f3-3*x2*x3*uy*f4+4*x2*x3*uy*f2+4*x2*ux*y3
279  & *f3+6*x2*ux*y3*f4-2*x2*ux*y3*f2-2*x2*ux*y2*f3+15*x2*
280  & ux*y2*f4-5*x2*ux*y2*f2+2*x3**2*uy*f3-6*x3**2*uy*f4+8*
281  & x3**2*uy*f2-2*x3*ux*y3*f3+6*x3*ux*y3*f4-8*x3*ux*y3*f2+
282  & x3*ux*y2*f3-3*x3*ux*y2*f4-2*x3*ux*y2*f2))*aux108
283  a23(ielem)=((8*x2**2*uy*f3-15*x2**2*uy*f4+5*x2**2*uy*f2
284  & -20*x2*x3*uy*f3+33*x2*x3*uy*f4-14*x2*x3*uy*f2+16*x2*
285  & ux*y3*f3-21*x2*ux*y3*f4+7*x2*ux*y3*f2-8*x2*ux*y2*f3+
286  & 15*x2*ux*y2*f4-5*x2*ux*y2*f2+8*x3**2*uy*f3-24*x3**2*uy
287  & *f4+17*x3**2*uy*f2-8*x3*ux*y3*f3+24*x3*ux*y3*f4-17*x3
288  & *ux*y3*f2+4*x3*ux*y2*f3-12*x3*ux*y2*f4+7*x3*ux*y2*f2))*aux108
289  a31(ielem)=((8*x2**2*uy*f3-6*x2**2*uy*f4+2*x2**2*uy*f2+
290  & 4*x2*x3*uy*f3-3*x2*x3*uy*f4-5*x2*x3*uy*f2-2*x2*ux*y3*
291  & f3-3*x2*ux*y3*f4+x2*ux*y3*f2-8*x2*ux*y2*f3+6*x2*ux*y2*
292  & f4-2*x2*ux*y2*f2+5*x3**2*uy*f3-15*x3**2*uy*f4+2*x3**2
293  & *uy*f2-5*x3*ux*y3*f3+15*x3*ux*y3*f4-2*x3*ux*y3*f2-2*
294  & x3*ux*y2*f3+6*x3*ux*y2*f4+4*x3*ux*y2*f2))*aux108
295  a32(ielem)=((17*x2**2*uy*f3-24*x2**2*uy*f4+8*x2**2*uy*
296  & f2-14*x2*x3*uy*f3+33*x2*x3*uy*f4-20*x2*x3*uy*f2+7*x2*
297  & ux*y3*f3-12*x2*ux*y3*f4+4*x2*ux*y3*f2-17*x2*ux*y2*f3+
298  & 24*x2*ux*y2*f4-8*x2*ux*y2*f2+5*x3**2*uy*f3-15*x3**2*uy
299  & *f4+8*x3**2*uy*f2-5*x3*ux*y3*f3+15*x3*ux*y3*f4-8*x3*
300  & ux*y3*f2+7*x3*ux*y2*f3-21*x3*ux*y2*f4+16*x3*ux*y2*f2))*aux108
301  a41(ielem)=(-(2*x2**2*uy*f3-15*x2**2*uy*f4+5*x2**2*uy*
302  & f2+x2*x3*uy*f3+6*x2*x3*uy*f4+x2*x3*uy*f2-2*x2*ux*y3*f3-
303  & 3*x2*ux*y3*f4+x2*ux*y3*f2-2*x2*ux*y2*f3+15*x2*ux*y2*f4-
304  & 5*x2*ux*y2*f2+5*x3**2*uy*f3-15*x3**2*uy*f4+2*x3**2*uy*
305  & f2-5*x3*ux*y3*f3+15*x3*ux*y3*f4-2*x3*ux*y3*f2+x3*ux*y2
306  & *f3-3*x3*ux*y2*f4-2*x3*ux*y2*f2))*aux036
307  a42(ielem)=(-(8*x2**2*uy*f3-24*x2**2*uy*f4+8*x2**2*uy*
308  & f2-11*x2*x3*uy*f3+24*x2*x3*uy*f4-8*x2*x3*uy*f2+7*x2*
309  & ux*y3*f3-12*x2*ux*y3*f4+4*x2*ux*y3*f2-8*x2*ux*y2*f3+
310  & 24*x2*ux*y2*f4-8*x2*ux*y2*f2+5*x3**2*uy*f3-15*x3**2*uy
311  & *f4+8*x3**2*uy*f2-5*x3*ux*y3*f3+15*x3*ux*y3*f4-8*x3*
312  & ux*y3*f2+4*x3*ux*y2*f3-12*x3*ux*y2*f4+4*x3*ux*y2*f2))*aux036
313  a43(ielem)=(-(8*x2**2*uy*f3-15*x2**2*uy*f4+5*x2**2*uy*
314  & f2-8*x2*x3*uy*f3+24*x2*x3*uy*f4-11*x2*x3*uy*f2+4*x2*
315  & ux*y3*f3-12*x2*ux*y3*f4+4*x2*ux*y3*f2-8*x2*ux*y2*f3+
316  & 15*x2*ux*y2*f4-5*x2*ux*y2*f2+8*x3**2*uy*f3-24*x3**2*uy
317  & *f4+8*x3**2*uy*f2-8*x3*ux*y3*f3+24*x3*ux*y3*f4-8*x3*
318  & ux*y3*f2+4*x3*ux*y2*f3-12*x3*ux*y2*f4+7*x3*ux*y2*f2))*aux036
319 !
320 ! DIAGONAL TERMS
321 ! THE SUM OF EACH COLUMN IS 0
322 !
323  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
324  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
325  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
326 !
327  ENDDO ! IELEM
328 !
329  ELSE
330 !
331  WRITE(lu,201) icoord
332  CALL plante(0)
333  stop
334  ENDIF
335 !
336 !
337  ELSEIF(ielmu.EQ.11) THEN
338 !
339 !================================
340 ! DERIVATIVE WRT X =
341 !================================
342 !
343  IF(icoord.EQ.1) THEN
344 !
345 ! LOOP ON THE ELEMENTS
346 !
347  DO ielem = 1 , nelem
348 !
349 ! INITIALISES THE GEOMETRICAL VARIABLES
350 !
351  x2 = xel(ielem,2)
352  x3 = xel(ielem,3)
353  y2 = yel(ielem,2)
354  y3 = yel(ielem,3)
355 !
356  f1 = f(ikle1(ielem))
357  f2 = f(ikle2(ielem)) - f1
358  f3 = f(ikle3(ielem)) - f1
359  f4 = f(ikle4(ielem)) - f1
360 !
361  u1 = u(ikle1(ielem))
362  u2 = u(ikle2(ielem))
363  u3 = u(ikle3(ielem))
364  v1 = v(ikle1(ielem))
365  v2 = v(ikle2(ielem))
366  v3 = v(ikle3(ielem))
367 !
368  ax1296 = xs1296 / surfac(ielem)
369  aux432 = xsu432 / surfac(ielem)
370 !
371 ! EXTRADIAGONAL TERMS
372 !
373  a12(ielem)=((((f3-3*f4)*y3+y2*f3)*(5*v3+2*v2+5*v1)*x2-
374  & 2*((f3-3*f4)*y3+y2*f3)*(5*v3+2*v2+5*v1)*x3+2*((3*f4
375  & -f2)*y2-y3*f2)*(5*v3+26*v2+17*v1)*x2-((3*f4-f2)*y2-y3
376  & *f2)*(5*v3+26*v2+17*v1)*x3+(y3*f3-3*y3*f4+y2*f3)*(2*
377  & y3-y2)*(5*u3+2*u2+5*u1)-(y3*f2-3*y2*f4+y2*f2)*(y3-2*
378  & y2)*(5*u3+26*u2+17*u1)))*ax1296
379  a13(ielem)=((((f3-3*f4)*y3+y2*f3)*(26*v3+5*v2+17*v1)*
380  & x2-2*((f3-3*f4)*y3+y2*f3)*(26*v3+5*v2+17*v1)*x3+2*(
381  & (3*f4-f2)*y2-y3*f2)*(2*v3+5*v2+5*v1)*x2-((3*f4-f2)*
382  & y2-y3*f2)*(2*v3+5*v2+5*v1)*x3+(y3*f3-3*y3*f4+y2*f3)*(
383  & 2*y3-y2)*(26*u3+5*u2+17*u1)-(y3*f2-3*y2*f4+y2*f2)*(y3
384  & -2*y2)*(2*u3+5*u2+5*u1)))*ax1296
385  a21(ielem)=(-(((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*(
386  & 5*v3+5*v2+2*v1)*x2-2*((2*f3-3*f4+f2)*y2-(f3-3*f4+2
387  & *f2)*y3)*(5*v3+5*v2+2*v1)*x3-((3*f4-f2)*y2-y3*f2)*(5
388  & *v3+17*v2+26*v1)*x2-((3*f4-f2)*y2-y3*f2)*(5*v3+17*v2
389  & +26*v1)*x3-(y3*f3-3*y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-y2
390  & *f2)*(2*y3-y2)*(5*u3+5*u2+2*u1)-(y3*f2-3*y2*f4+y2*f2
391  & )*(y3+y2)*(5*u3+17*u2+26*u1)))*ax1296
392  a23(ielem)=(-(((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*(
393  & 26*v3+17*v2+5*v1)*x2-2*((2*f3-3*f4+f2)*y2-(f3-3*f4+
394  & 2*f2)*y3)*(26*v3+17*v2+5*v1)*x3-((3*f4-f2)*y2-y3*f2)*
395  & (2*v3+5*v2+5*v1)*x2-((3*f4-f2)*y2-y3*f2)*(2*v3+5*v2
396  & +5*v1)*x3-(y3*f3-3*y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-y2*
397  & f2)*(2*y3-y2)*(26*u3+17*u2+5*u1)-(y3*f2-3*y2*f4+y2*
398  & f2)*(y3+y2)*(2*u3+5*u2+5*u1)))*ax1296
399  a31(ielem)=(-(2*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)
400  & *(5*v3+5*v2+2*v1)*x2-((2*f3-3*f4+f2)*y2-(f3-3*f4+2
401  & *f2)*y3)*(5*v3+5*v2+2*v1)*x3+((f3-3*f4)*y3+y2*f3)*(
402  & 17*v3+5*v2+26*v1)*x2+((f3-3*f4)*y3+y2*f3)*(17*v3+5*
403  & v2+26*v1)*x3-(y3*f3-3*y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-
404  & y2*f2)*(y3-2*y2)*(5*u3+5*u2+2*u1)-(y3*f3-3*y3*f4+y2*
405  & f3)*(y3+y2)*(17*u3+5*u2+26*u1)))*ax1296
406  a32(ielem)=(-(2*((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)
407  & *(17*v3+26*v2+5*v1)*x2-((2*f3-3*f4+f2)*y2-(f3-3*f4+
408  & 2*f2)*y3)*(17*v3+26*v2+5*v1)*x3+((f3-3*f4)*y3+y2*f3)*
409  & (5*v3+2*v2+5*v1)*x2+((f3-3*f4)*y3+y2*f3)*(5*v3+2*v2
410  & +5*v1)*x3-(y3*f3-3*y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-y2*
411  & f2)*(y3-2*y2)*(17*u3+26*u2+5*u1)-(y3*f3-3*y3*f4+y2*
412  & f3)*(y3+y2)*(5*u3+2*u2+5*u1)))*ax1296
413  a41(ielem)=((((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*(5
414  & *v3+5*v2+2*v1)*x2-((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)
415  & *y3)*(5*v3+5*v2+2*v1)*x3+((f3-3*f4)*y3+y2*f3)*(17*v3
416  & +5*v2+26*v1)*x3-((f3-3*f4)*y3+y2*f3)*(17*u3+5*u2+26
417  & *u1)*y3-((3*f4-f2)*y2-y3*f2)*(5*v3+17*v2+26*v1)*x2+((
418  & 3*f4-f2)*y2-y3*f2)*(5*u3+17*u2+26*u1)*y2-(y3*f3-3*y3*
419  & f4+2*y3*f2-2*y2*f3+3*y2*f4-y2*f2)*(y3-y2)*(5*u3+5*u2
420  & +2*u1)))*aux432
421  a42(ielem)=((((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*(
422  & 17*v3+26*v2+5*v1)*x2-((2*f3-3*f4+f2)*y2-(f3-3*f4+2*
423  & f2)*y3)*(17*v3+26*v2+5*v1)*x3+((f3-3*f4)*y3+y2*f3)*(
424  & 5*v3+2*v2+5*v1)*x3-((f3-3*f4)*y3+y2*f3)*(5*u3+2*u2+
425  & 5*u1)*y3-((3*f4-f2)*y2-y3*f2)*(5*v3+26*v2+17*v1)*x2+(
426  & (3*f4-f2)*y2-y3*f2)*(5*u3+26*u2+17*u1)*y2-(y3*f3-3*
427  & y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-y2*f2)*(y3-y2)*(17*u3+
428  & 26*u2+5*u1)))*aux432
429  a43(ielem)=((((2*f3-3*f4+f2)*y2-(f3-3*f4+2*f2)*y3)*(
430  & 26*v3+17*v2+5*v1)*x2-((2*f3-3*f4+f2)*y2-(f3-3*f4+2*
431  & f2)*y3)*(26*v3+17*v2+5*v1)*x3+((f3-3*f4)*y3+y2*f3)*(
432  & 26*v3+5*v2+17*v1)*x3-((f3-3*f4)*y3+y2*f3)*(26*u3+5*
433  & u2+17*u1)*y3-((3*f4-f2)*y2-y3*f2)*(2*v3+5*v2+5*v1)*
434  & x2+((3*f4-f2)*y2-y3*f2)*(2*u3+5*u2+5*u1)*y2-(y3*f3-3
435  & *y3*f4+2*y3*f2-2*y2*f3+3*y2*f4-y2*f2)*(y3-y2)*(26*u3+
436  & 17*u2+5*u1)))*aux432
437 !
438 ! DIAGONAL TERMS
439 ! THE SUM OF EACH COLUMN IS 0
440 !
441  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
442  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
443  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
444 !
445  ENDDO ! IELEM
446 !
447  ELSEIF(icoord.EQ.2) THEN
448 !
449 !================================
450 ! DERIVATIVE WRT Y =
451 !================================
452 !
453  DO ielem = 1 , nelem
454 !
455 ! INITIALISES THE GEOMETRICAL VARIABLES
456 !
457  x2 = xel(ielem,2)
458  x3 = xel(ielem,3)
459  y2 = yel(ielem,2)
460  y3 = yel(ielem,3)
461 !
462  f1 = f(ikle1(ielem))
463  f2 = f(ikle2(ielem)) - f1
464  f3 = f(ikle3(ielem)) - f1
465  f4 = f(ikle4(ielem)) - f1
466 !
467  u1 = u(ikle1(ielem))
468  u2 = u(ikle2(ielem))
469  u3 = u(ikle3(ielem))
470  v1 = v(ikle1(ielem))
471  v2 = v(ikle2(ielem))
472  v3 = v(ikle3(ielem))
473 !
474  ax1296 = xs1296 / surfac(ielem)
475  aux432 = xsu432 / surfac(ielem)
476 !
477 ! EXTRADIAGONAL TERMS
478 !
479  a12(ielem)=(-(((2*y3-y2)*(5*u3+2*u2+5*u1)*f3-(f3+3*f4
480  & )*(5*v3+2*v2+5*v1)*x3)*x2+((y3-2*y2)*(3*f4-f2)*(5*
481  & u3+26*u2+17*u1)-(3*f4+f2)*(5*v3+26*v2+17*v1)*x3)*x2
482  & +(2*y3-y2)*(f3-3*f4)*(5*u3+2*u2+5*u1)*x3-(y3-2*y2)*
483  & (5*u3+26*u2+17*u1)*x3*f2-2*(f3-3*f4)*(5*v3+2*v2+5
484  & *v1)*x3**2+2*(3*f4-f2)*(5*v3+26*v2+17*v1)*x2**2+(5*
485  & v3+26*v2+17*v1)*x3**2*f2+(5*v3+2*v2+5*v1)*x2**2*f3))*ax1296
486  a13(ielem)=(-(((2*y3-y2)*(26*u3+5*u2+17*u1)*f3-(f3+3*
487  & f4)*(26*v3+5*v2+17*v1)*x3)*x2+((y3-2*y2)*(3*f4-f2)*(
488  & 2*u3+5*u2+5*u1)-(3*f4+f2)*(2*v3+5*v2+5*v1)*x3)*x2+(
489  & 2*y3-y2)*(f3-3*f4)*(26*u3+5*u2+17*u1)*x3-(y3-2*y2)*(
490  & 2*u3+5*u2+5*u1)*x3*f2-2*(f3-3*f4)*(26*v3+5*v2+17*
491  & v1)*x3**2+2*(3*f4-f2)*(2*v3+5*v2+5*v1)*x2**2+(26*v3
492  & +5*v2+17*v1)*x2**2*f3+(2*v3+5*v2+5*v1)*x3**2*f2))*ax1296
493  a21(ielem)=((((2*y3-y2)*(2*f3-3*f4+f2)*(5*u3+5*u2+2*
494  & u1)-(5*f3-9*f4+4*f2)*(5*v3+5*v2+2*v1)*x3)*x2+((y3+
495  & y2)*(3*f4-f2)*(5*u3+17*u2+26*u1)-(3*f4-2*f2)*(5*v3
496  & +17*v2+26*v1)*x3)*x2-(2*y3-y2)*(f3-3*f4+2*f2)*(5*u3
497  & +5*u2+2*u1)*x3-(y3+y2)*(5*u3+17*u2+26*u1)*x3*f2+(2*
498  & f3-3*f4+f2)*(5*v3+5*v2+2*v1)*x2**2+2*(f3-3*f4+2*f2
499  & )*(5*v3+5*v2+2*v1)*x3**2-(3*f4-f2)*(5*v3+17*v2+26*
500  & v1)*x2**2+(5*v3+17*v2+26*v1)*x3**2*f2))*ax1296
501  a23(ielem)=((((2*y3-y2)*(2*f3-3*f4+f2)*(26*u3+17*u2+
502  & 5*u1)-(5*f3-9*f4+4*f2)*(26*v3+17*v2+5*v1)*x3)*x2+((
503  & y3+y2)*(3*f4-f2)*(2*u3+5*u2+5*u1)-(3*f4-2*f2)*(2*
504  & v3+5*v2+5*v1)*x3)*x2-(2*y3-y2)*(f3-3*f4+2*f2)*(26*
505  & u3+17*u2+5*u1)*x3-(y3+y2)*(2*u3+5*u2+5*u1)*x3*f2+(2
506  & *f3-3*f4+f2)*(26*v3+17*v2+5*v1)*x2**2+2*(f3-3*f4+2
507  & *f2)*(26*v3+17*v2+5*v1)*x3**2-(3*f4-f2)*(2*v3+5*v2+
508  & 5*v1)*x2**2+(2*v3+5*v2+5*v1)*x3**2*f2))*ax1296
509  a31(ielem)=(-(((y3+y2)*(17*u3+5*u2+26*u1)*f3-(2*f3-3*
510  & f4)*(17*v3+5*v2+26*v1)*x3)*x2-((y3-2*y2)*(2*f3-3*f4
511  & +f2)*(5*u3+5*u2+2*u1)-(4*f3-9*f4+5*f2)*(5*v3+5*v2
512  & +2*v1)*x3)*x2+(y3+y2)*(f3-3*f4)*(17*u3+5*u2+26*u1)*
513  & x3+(y3-2*y2)*(f3-3*f4+2*f2)*(5*u3+5*u2+2*u1)*x3-2*
514  & (2*f3-3*f4+f2)*(5*v3+5*v2+2*v1)*x2**2-(f3-3*f4+2*
515  & f2)*(5*v3+5*v2+2*v1)*x3**2-(f3-3*f4)*(17*v3+5*v2+
516  & 26*v1)*x3**2-(17*v3+5*v2+26*v1)*x2**2*f3))*ax1296
517  a32(ielem)=(-(((y3+y2)*(5*u3+2*u2+5*u1)*f3-(2*f3-3*f4
518  & )*(5*v3+2*v2+5*v1)*x3)*x2-((y3-2*y2)*(2*f3-3*f4+f2)
519  & *(17*u3+26*u2+5*u1)-(4*f3-9*f4+5*f2)*(17*v3+26*v2
520  & +5*v1)*x3)*x2+(y3+y2)*(f3-3*f4)*(5*u3+2*u2+5*u1)*x3+
521  & (y3-2*y2)*(f3-3*f4+2*f2)*(17*u3+26*u2+5*u1)*x3-2*(
522  & 2*f3-3*f4+f2)*(17*v3+26*v2+5*v1)*x2**2-(f3-3*f4+2*
523  & f2)*(17*v3+26*v2+5*v1)*x3**2-(f3-3*f4)*(5*v3+2*v2+
524  & 5*v1)*x3**2-(5*v3+2*v2+5*v1)*x2**2*f3))*ax1296
525  a41(ielem)=(-(((y3-y2)*(2*f3-3*f4+f2)*(5*u3+5*u2+2*u1
526  & )-3*(f3-2*f4+f2)*(5*v3+5*v2+2*v1)*x3)*x2+((3*f4-f2)
527  & *(5*u3+17*u2+26*u1)*y2+(5*v3+17*v2+26*v1)*x3*f2)*x2
528  & +((17*v3+5*v2+26*v1)*x3-(17*u3+5*u2+26*u1)*y3)*x2*
529  & f3-(y3-y2)*(f3-3*f4+2*f2)*(5*u3+5*u2+2*u1)*x3+(2*f3
530  & -3*f4+f2)*(5*v3+5*v2+2*v1)*x2**2+(f3-3*f4+2*f2)*(5
531  & *v3+5*v2+2*v1)*x3**2+(f3-3*f4)*(17*v3+5*v2+26*v1)*
532  & x3**2-(f3-3*f4)*(17*u3+5*u2+26*u1)*x3*y3-(3*f4-f2)*(
533  & 5*v3+17*v2+26*v1)*x2**2-(5*u3+17*u2+26*u1)*x3*y2*f2))*aux432
534  a42(ielem)=(-(((y3-y2)*(2*f3-3*f4+f2)*(17*u3+26*u2+5*
535  & u1)-3*(f3-2*f4+f2)*(17*v3+26*v2+5*v1)*x3)*x2+((3*f4
536  & -f2)*(5*u3+26*u2+17*u1)*y2+(5*v3+26*v2+17*v1)*x3*f2
537  & )*x2+((5*v3+2*v2+5*v1)*x3-(5*u3+2*u2+5*u1)*y3)*x2*
538  & f3-(y3-y2)*(f3-3*f4+2*f2)*(17*u3+26*u2+5*u1)*x3+(2*
539  & f3-3*f4+f2)*(17*v3+26*v2+5*v1)*x2**2+(f3-3*f4+2*f2)
540  & *(17*v3+26*v2+5*v1)*x3**2+(f3-3*f4)*(5*v3+2*v2+5*
541  & v1)*x3**2-(f3-3*f4)*(5*u3+2*u2+5*u1)*x3*y3-(3*f4-f2)
542  & *(5*v3+26*v2+17*v1)*x2**2-(5*u3+26*u2+17*u1)*x3*y2*
543  & f2))*aux432
544  a43(ielem)=(-(((y3-y2)*(2*f3-3*f4+f2)*(26*u3+17*u2+5*
545  & u1)-3*(f3-2*f4+f2)*(26*v3+17*v2+5*v1)*x3)*x2+((3*f4
546  & -f2)*(2*u3+5*u2+5*u1)*y2+(2*v3+5*v2+5*v1)*x3*f2)*x2
547  & +((26*v3+5*v2+17*v1)*x3-(26*u3+5*u2+17*u1)*y3)*x2*
548  & f3-(y3-y2)*(f3-3*f4+2*f2)*(26*u3+17*u2+5*u1)*x3+(2*
549  & f3-3*f4+f2)*(26*v3+17*v2+5*v1)*x2**2+(f3-3*f4+2*f2)
550  & *(26*v3+17*v2+5*v1)*x3**2+(f3-3*f4)*(26*v3+5*v2+17
551  & *v1)*x3**2-(f3-3*f4)*(26*u3+5*u2+17*u1)*x3*y3-(3*f4-
552  & f2)*(2*v3+5*v2+5*v1)*x2**2-(2*u3+5*u2+5*u1)*x3*y2*f2))*aux432
553 !
554 ! DIAGONAL TERMS
555 ! THE SUM OF EACH COLUMN IS 0
556 !
557  a11(ielem) = - a21(ielem) - a31(ielem) - a41(ielem)
558  a22(ielem) = - a12(ielem) - a32(ielem) - a42(ielem)
559  a33(ielem) = - a13(ielem) - a23(ielem) - a43(ielem)
560 !
561  ENDDO ! IELEM
562 !
563  ELSE
564 !
565  WRITE(lu,201) icoord
566  CALL plante(0)
567  stop
568  ENDIF
569 !
570  ELSE
571 !
572  WRITE(lu,301) ielmu
573  CALL plante(0)
574  stop
575 !
576  ENDIF
577 !
578 !-----------------------------------------------------------------------
579 !
580  ELSE
581  WRITE(lu,101) ielmf
582 101 FORMAT(1x,'MT12BA (BIEF) :',/,
583  & 1x,'DISCRETIZATION OF F : ',1i6,' NOT AVAILABLE')
584  CALL plante(0)
585  stop
586  ENDIF
587 !
588 201 FORMAT(1x,'MT12BA (BIEF) : IMPOSSIBLE COMPONENT ',
589  & 1i6,' CHECK ICOORD')
590 !
591 301 FORMAT(1x,'MT12BA (BIEF) :',/,
592  & 1x,'DISCRETIZATION OF U : ',1i6,' NOT AVAILABLE')
593 !
594 !-----------------------------------------------------------------------
595 !
596  RETURN
597  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mt12ba
subroutine mt12ba(A11, A12, A13, A21, A22, A23, A31, A32, A33, A41, A42, A43, XMUL, SF, SU, SV, F, U, V, XEL, YEL, SURFAC, IKLE1, IKLE2, IKLE3, IKLE4, NELEM, NELMAX, ICOORD)
Definition: mt12ba.f:14
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:13 for The TELEMAC-MASCARET system by doxygen 1.8.13