The TELEMAC-MASCARET system  trunk
mt02bb.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt02bb
3 ! *****************
4 !
5  &( a11 , a12 , a13 , a14 ,
6  & a22 , a23 , a24 ,
7  & a33 , a34 ,
8  & a44 ,
9  & xmul,su,u,xel,yel,surfac,ikle1,ikle2,ikle3,nelem,nelmax)
10 !
11 !***********************************************************************
12 ! BIEF V6P1 21/08/2010
13 !***********************************************************************
14 !
15 !brief BUILDS THE DIFFUSION MATRIX FOR QUASI-BUBBLE TRIANGLES.
16 !+
17 !+ VISCOSITY CAN BE ISOTROPIC, OR NOT ISOTROPIC. IN THIS
18 !+ CASE U IS AN ARRAY WITH SECOND DIMENSION EQUAL TO 3.
19 !
20 !warning THE JACOBIAN MUST BE POSITIVE
21 !
22 !history J-M HERVOUET (LNH) ; C MOULIN (LNH)
23 !+ 10/01/95
24 !+ V5P1
25 !+
26 !
27 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
28 !+ 13/07/2010
29 !+ V6P0
30 !+ Translation of French comments within the FORTRAN sources into
31 !+ English comments
32 !
33 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
34 !+ 21/08/2010
35 !+ V6P0
36 !+ Creation of DOXYGEN tags for automated documentation and
37 !+ cross-referencing of the FORTRAN sources
38 !
39 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
40 !| A11 |<--| ELEMENTS OF MATRIX
41 !| A12 |<--| ELEMENTS OF MATRIX
42 !| A13 |<--| ELEMENTS OF MATRIX
43 !| A14 |<--| ELEMENTS OF MATRIX
44 !| A22 |<--| ELEMENTS OF MATRIX
45 !| A23 |<--| ELEMENTS OF MATRIX
46 !| A24 |<--| ELEMENTS OF MATRIX
47 !| A33 |<--| ELEMENTS OF MATRIX
48 !| A34 |<--| ELEMENTS OF MATRIX
49 !| A44 |<--| ELEMENTS OF MATRIX
50 !| IKLE1 |-->| FIRST POINTS OF TRIANGLES
51 !| IKLE2 |-->| SECOND POINTS OF TRIANGLES
52 !| IKLE3 |-->| THIRD POINTS OF TRIANGLES
53 !| NELEM |-->| NUMBER OF ELEMENTS
54 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
55 !| SU |-->| BIEF_OBJ STRUCTURE OF U
56 !| SURFAC |-->| AREA OF TRIANGLES
57 !| U |-->| FUNCTION U USED IN THE FORMULA
58 !| XEL |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
59 !| YEL |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
60 !| XMUL |-->| MULTIPLICATION FACTOR
61 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
62 !
63  USE bief, ex_mt02bb => mt02bb
64 !
65  USE declarations_special
66  IMPLICIT NONE
67 !
68 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
69 !
70  INTEGER, INTENT(IN) :: NELEM,NELMAX
71  INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax),IKLE3(nelmax)
72  DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*),A14(*)
73  DOUBLE PRECISION, INTENT(INOUT) :: A22(*),A23(*),A24(*)
74  DOUBLE PRECISION, INTENT(INOUT) :: A33(*),A34(*)
75  DOUBLE PRECISION, INTENT(INOUT) :: A44(*)
76  DOUBLE PRECISION, INTENT(IN) :: XMUL,U(*)
77 ! STRUCTURE OF U
78  TYPE(bief_obj), INTENT(IN) :: SU
79 !
80  DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
81  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
82 !
83 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
84 !
85 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
86 !
87  INTEGER IELMNU,IELEM,ISO,IAD2,IAD3
88 !
89  DOUBLE PRECISION X2,X3,Y2,Y3,AUX1,AUX2
90  DOUBLE PRECISION NUX1,NUX2,NUX3,NUY1,NUY2,NUY3,NUZ1,NUZ2,NUZ3
91 !
92 !=======================================================================
93 !
94 ! EXTRACTS THE TYPE OF ELEMENT FOR VISCOSITY
95 !
96  ielmnu = su%ELM
97  iso = su%DIM2
98 !
99 ! IF(IELMNU.EQ.10.AND.ISO.EQ.1) THEN
100 !
101 !-----------------------------------------------------------------------
102 !
103 ! P0 DISCRETISATION FOR VISCOSITY:
104 !
105 !-----------------------------------------------------------------------
106 !
107  IF(ielmnu.EQ.11.AND.iso.EQ.1) THEN
108 !
109 !-----------------------------------------------------------------------
110 !
111 ! P1 DISCRETISATION FOR ISOTROPIC VISCOSITY:
112 !
113  DO ielem = 1 , nelem
114 !
115 ! INITIALISES THE GEOMETRICAL VARIABLES
116 !
117  x2 = xel(ielem,2)
118  x3 = xel(ielem,3)
119 !
120  y2 = yel(ielem,2)
121  y3 = yel(ielem,3)
122 !
123  nux1 = u(ikle1(ielem))
124  nux2 = u(ikle2(ielem))
125  nux3 = u(ikle3(ielem))
126 !
127  aux1 = xmul/(108*surfac(ielem))
128  aux2 = 3 * aux1
129 !
130 ! EXTRADIAGONAL TERMS
131 ! SOME FACTORISATIONS REMAIN TO BE DONE
132 !
133  a12(ielem)=((2*x2**2+x2*x3-x3**2-y3**2+y3*y2+2*y2**2)*(
134  & nux3+4*nux2+4*nux1))*aux1
135 !
136  a13(ielem)=(-(x2**2-x2*x3-2*x3**2-2*y3**2-y3*y2+y2**2)*(
137  & 4*nux3+nux2+4*nux1))*aux1
138 !
139  a14(ielem)=((4*nux3+nux2+4*nux1)*x2*x3-2*(4*nux3+nux2+
140  & 4*nux1)*x3**2-2*(4*nux3+nux2+4*nux1)*y3**2+(4*nux3+
141  & nux2+4*nux1)*y3*y2-2*(nux3+4*nux2+4*nux1)*x2**2+(nux3
142  & +4*nux2+4*nux1)*x2*x3+(nux3+4*nux2+4*nux1)*y3*y2-2*(
143  & nux3+4*nux2+4*nux1)*y2**2)*aux2
144 !
145  a23(ielem)=((2*x2**2-5*x2*x3+2*x3**2+2*y3**2-5*y3*y2+
146  & 2*y2**2)*(4*nux3+4*nux2+nux1))*aux1
147 !
148  a24(ielem)=(-((4*nux3+4*nux2+nux1)*x2**2-3*(4*nux3+4*
149  & nux2+nux1)*x2*x3+2*(4*nux3+4*nux2+nux1)*x3**2+2*(4*
150  & nux3+4*nux2+nux1)*y3**2-3*(4*nux3+4*nux2+nux1)*y3*y2+
151  & (4*nux3+4*nux2+nux1)*y2**2+(nux3+4*nux2+4*nux1)*x2**2
152  & +(nux3+4*nux2+4*nux1)*x2*x3+(nux3+4*nux2+4*nux1)*y3*
153  & y2+(nux3+4*nux2+4*nux1)*y2**2))*aux2
154 !
155  a34(ielem)=(
156  & nux1*(-5*y3**2-y3*y2-2*y2**2)+nux2*(-5*y3**
157  & 2+11*y3*y2-8*y2**2)+8*nux3*(-y3**2+y3*y2-y2**2)+nux1*(
158  & -2*x2**2-x2*x3-5*x3**2)+nux2*(-8*x2**2+11*x2*x3-5*x3
159  & **2)+8*nux3*(-x2**2+x2*x3-x3**2))*aux2
160 !
161 ! THE DIAGONAL TERMS ARE OBTAINED BY MEANS OF THE
162 ! MAGIC SQUARE:
163 !
164  a11(ielem) = - a12(ielem) - a13(ielem) - a14(ielem)
165  a22(ielem) = - a12(ielem) - a23(ielem) - a24(ielem)
166  a33(ielem) = - a13(ielem) - a23(ielem) - a34(ielem)
167  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
168 !
169  ENDDO ! IELEM
170 !
171 !-----------------------------------------------------------------------
172 !
173  ELSEIF(ielmnu.EQ.11.AND.iso.EQ.3) THEN
174 !
175 !-----------------------------------------------------------------------
176 !
177 ! P1 DISCRETISATION FOR NONISOTROPIC VISCOSITY:
178 !
179  iad2 = su%MAXDIM1
180  iad3 = 2*iad2
181  DO ielem = 1 , nelem
182 !
183 ! INITIALISES THE GEOMETRICAL VARIABLES
184 !
185  x2 = xel(ielem,2)
186  x3 = xel(ielem,3)
187 !
188  y2 = yel(ielem,2)
189  y3 = yel(ielem,3)
190 !
191  nux1 = u(ikle1(ielem))
192  nux2 = u(ikle2(ielem))
193  nux3 = u(ikle3(ielem))
194  nuy1 = u(ikle1(ielem)+iad2)
195  nuy2 = u(ikle2(ielem)+iad2)
196  nuy3 = u(ikle3(ielem)+iad2)
197  nuz1 = u(ikle1(ielem)+iad3)
198  nuz2 = u(ikle2(ielem)+iad3)
199  nuz3 = u(ikle3(ielem)+iad3)
200 !
201  aux1 = xmul/(108*surfac(ielem))
202  aux2 = 3 * aux1
203 !
204 ! EXTRADIAGONAL TERMS
205 !
206  a12(ielem)=(-4*nux1*(y3+y2)*(y3-2*y2)-4*nux2*(y3+y2)*(
207  & y3-2*y2)-nux3*(y3+y2)*(y3-2*y2)+((nuy3+4*nuy2+4*nuy1)
208  & *x3-(y3+4*y2)*nuz3-4*(y3+4*y2)*nuz2-4*(y3+4*y2)*nuz1
209  & )*x2+(nuz3+4*nuz2+4*nuz1)*(2*y3-y2)*x3+2*(nuy3+4*
210  & nuy2+4*nuy1)*x2**2-(nuy3+4*nuy2+4*nuy1)*x3**2)*aux1
211 !
212  a13(ielem)=(4*nux1*(2*y3-y2)*(y3+y2)+nux2*(2*y3-y2)*(y3
213  & +y2)+4*nux3*(2*y3-y2)*(y3+y2)+((4*nuy3+nuy2+4*nuy1)*
214  & x3-4*(y3-2*y2)*nuz3-(y3-2*y2)*nuz2-4*(y3-2*y2)*nuz1)
215  & *x2-(4*nuz3+nuz2+4*nuz1)*(4*y3+y2)*x3-(4*nuy3+nuy2+4
216  & *nuy1)*x2**2+2*(4*nuy3+nuy2+4*nuy1)*x3**2)*aux1
217 !
218  a14(ielem)=(4*nux1*(-(2*y3-y2)*y3+(y3-2*y2)*y2)+nux2*(-
219  & (2*y3-y2)*y3+4*(y3-2*y2)*y2)+nux3*(-4*(2*y3-y2)*y3+(
220  & y3-2*y2)*y2)+((4*nuy3+nuy2+4*nuy1)*x3-4*nuz3*y3-nuz2*
221  & y3-4*nuz1*y3)*x2+((nuy3+4*nuy2+4*nuy1)*x3-(y3-4*y2)*
222  & nuz3-4*(y3-4*y2)*nuz2-4*(y3-4*y2)*nuz1)*x2+(4*nuz3+
223  & nuz2+4*nuz1)*(4*y3-y2)*x3-(nuz3+4*nuz2+4*nuz1)*x3*y2-
224  & 2*(4*nuy3+nuy2+4*nuy1)*x3**2-2*(nuy3+4*nuy2+4*nuy1)*x2**2)*aux2
225 !
226  a23(ielem)=(-((5*(4*nuy3+4*nuy2+nuy1)*x3-4*(5*y3-4*
227  & y2)*nuz3-4*(5*y3-4*y2)*nuz2-(5*y3-4*y2)*nuz1)*x2+(4
228  & *nuz3+4*nuz2+nuz1)*(4*y3-5*y2)*x3-2*(4*nuy3+4*nuy2+
229  & nuy1)*x2**2-2*(4*nuy3+4*nuy2+nuy1)*x3**2-4*(2*y3-y2)
230  & *(y3-2*y2)*nux3-4*(2*y3-y2)*(y3-2*y2)*nux2-(2*y3-y2)
231  & *(y3-2*y2)*nux1))*aux1
232 !
233  a24(ielem)=(-(5*x2**2*nuy3+8*x2**2*nuy2+5*x2**2*nuy1-
234  & 11*x2*x3*nuy3-8*x2*x3*nuy2+x2*x3*nuy1+11*x2*nuz3*y3-10
235  & *x2*nuz3*y2+8*x2*nuz2*y3-16*x2*nuz2*y2-x2*nuz1*y3-10*
236  & x2*nuz1*y2+8*x3**2*nuy3+8*x3**2*nuy2+2*x3**2*nuy1-16*
237  & x3*nuz3*y3+11*x3*nuz3*y2-16*x3*nuz2*y3+8*x3*nuz2*y2-4
238  & *x3*nuz1*y3-x3*nuz1*y2+8*nux3*y3**2-11*nux3*y3*y2+5*
239  & nux3*y2**2+8*nux2*y3**2-8*nux2*y3*y2+8*nux2*y2**2+2*
240  & nux1*y3**2+nux1*y3*y2+5*nux1*y2**2))*aux2
241 !
242  a34(ielem)=(-(8*x2**2*nuy3+8*x2**2*nuy2+2*x2**2*nuy1-8
243  & *x2*x3*nuy3-11*x2*x3*nuy2+x2*x3*nuy1+8*x2*nuz3*y3-16*
244  & x2*nuz3*y2+11*x2*nuz2*y3-16*x2*nuz2*y2-x2*nuz1*y3-4*x2
245  & *nuz1*y2+8*x3**2*nuy3+5*x3**2*nuy2+5*x3**2*nuy1-16*x3
246  & *nuz3*y3+8*x3*nuz3*y2-10*x3*nuz2*y3+11*x3*nuz2*y2-10*
247  & x3*nuz1*y3-x3*nuz1*y2+8*nux3*y3**2-8*nux3*y3*y2+8*nux3
248  & *y2**2+5*nux2*y3**2-11*nux2*y3*y2+8*nux2*y2**2+5*nux1
249  & *y3**2+nux1*y3*y2+2*nux1*y2**2))*aux2
250 !
251 ! THE DIAGONAL TERMS ARE OBTAINED BY MEANS OF THE
252 ! MAGIC SQUARE:
253 !
254  a11(ielem) = - a12(ielem) - a13(ielem) - a14(ielem)
255  a22(ielem) = - a12(ielem) - a23(ielem) - a24(ielem)
256  a33(ielem) = - a13(ielem) - a23(ielem) - a34(ielem)
257  a44(ielem) = - a14(ielem) - a24(ielem) - a34(ielem)
258 !
259  ENDDO ! IELEM
260 !
261 !-----------------------------------------------------------------------
262 !
263  ELSE
264 !
265  WRITE(lu,11) ielmnu,iso
266 11 FORMAT(1x,
267  & 'MT02BB (BIEF) : TYPE OF VISCOSITY NOT AVAILABLE : ',2i6)
268  CALL plante(0)
269  stop
270 !
271  ENDIF
272 !
273 !-----------------------------------------------------------------------
274 !
275  RETURN
276  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mt02bb
subroutine mt02bb(A11, A12, A13, A14, A22, A23, A24, A33, A34, A44, XMUL, SU, U, XEL, YEL, SURFAC, IKLE1, IKLE2, IKLE3, NELEM, NELMAX)
Definition: mt02bb.f:11
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:36:45 for The TELEMAC-MASCARET system by doxygen 1.8.13