The TELEMAC-MASCARET system  trunk
mt04pp.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt04pp
3 ! *****************
4 !
5  &( t,xm,xmul,su,sv,sw,u,v,x,y,z,surfac,ikle,nelem,nelmax,formul)
6 !
7 !***********************************************************************
8 ! BIEF V6P3 21/08/2010
9 !***********************************************************************
10 !
11 !brief BUILDS THE MATRIX U GRAG (PSII) U GRAD (PSIJ).
12 !code
13 !+ COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
14 !+
15 !+ MAUGUG2 :
16 !+
17 !+ / -> ->
18 !+ A = XMUL / U . GRAD (P ) * U . GRAD (P ) * J(X,Y) DXDY
19 !+ I J /S 2D 2D I 2D 2D J
20 !+
21 !+
22 !+ MAUGUG1 : SEE COMPONENTS F0 AND G0
23 !+
24 !+ / -> ->
25 !+ A = XMUL / F . GRAD (P ) * U . GRAD (P ) * J(X,Y) DXDY
26 !+ I J /S 2D 2D I 2D 2D J
27 !+
28 !+
29 !+ BY ELEMENTARY CELL; THE ELEMENT IS THE P1 PRISM
30 !+
31 !+ J(X,Y): JACOBIAN OF THE ISOPARAMETRIC TRANSFORMATION
32 !
33 !note SIMPLIFICATIONS: AN AVERAGE HEIGHT OF THE PRISM IS TAKEN.
34 !+ VELOCITIES ARE CONSIDERED CONSTANT ON THE ELEMENT.
35 !
36 !warning THE JACOBIAN MUST BE POSITIVE
37 !
38 !history J-M HERVOUET (LNHE) ; F LEPEINTRE (LNH)
39 !+ 07/07/2005
40 !+
41 !+
42 !
43 !history
44 !+ 16/09/2005
45 !+ V5P6
46 !+ VERTICAL UPWIND REMOVED
47 !
48 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
49 !+ 13/07/2010
50 !+ V6P0
51 !+ Translation of French comments within the FORTRAN sources into
52 !+ English comments
53 !
54 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
55 !+ 21/08/2010
56 !+ V6P0
57 !+ Creation of DOXYGEN tags for automated documentation and
58 !+ cross-referencing of the FORTRAN sources
59 !
60 !history J-M HERVOUET (EDF R&D, LNHE)
61 !+ 11/01/2013
62 !+ V6P3
63 !+ XEL and YEL sent instead of XPT and YPT for X and Y.
64 !
65 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
66 !| FORMUL |-->| FORMULA DESCRIBING THE RESULTING MATRIX
67 !| IKLE |-->| CONNECTIVITY TABLE.
68 !| NELEM |-->| NUMBER OF ELEMENTS
69 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
70 !| SU |-->| STRUCTURE OF FUNCTIONS U
71 !| SV |-->| STRUCTURE OF FUNCTIONS V
72 !| SW |-->| STRUCTURE OF FUNCTIONS W
73 !| SURFAC |-->| AREA OF 2D ELEMENTS
74 !| T |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
75 !| U |-->| FUNCTION USED IN THE FORMULA
76 !| V |-->| FUNCTION USED IN THE FORMULA
77 !| X |-->| ABSCISSAE OF POINTS, PER ELEMENT
78 !| Y |-->| ORDINATES OF POINTS, PER ELEMENT
79 !| Z |-->| ELEVATIONS OF POINTS, STILL PER POINTS !!!
80 !| XM |<->| OFF-DIAGONAL TERMS
81 !| XMUL |-->| COEFFICIENT FOR MULTIPLICATION
82 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
83 !
84  USE bief, ex_mt04pp => mt04pp
85 !
86  USE declarations_special
87  IMPLICIT NONE
88 !
89 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
90 !
91  INTEGER, INTENT(IN) :: NELEM,NELMAX
92  INTEGER, INTENT(IN) :: IKLE(nelmax,6)
93 !
94  DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,6),XM(nelmax,30)
95  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
96 !
97  DOUBLE PRECISION, INTENT(IN) :: XMUL
98  DOUBLE PRECISION, INTENT(IN) :: U(*),V(*)
99 !
100 ! STRUCTURES OF U, V, W
101  TYPE(bief_obj), INTENT(IN) :: SU,SV,SW
102 !
103  DOUBLE PRECISION, INTENT(IN) :: X(nelmax,6),Y(nelmax,6),Z(*)
104 !
105  CHARACTER(LEN=16), INTENT(IN) :: FORMUL
106 !
107 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
108 !
109 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
110 !
111  DOUBLE PRECISION X2,X3,Y2,Y3
112  DOUBLE PRECISION U0,V0,F0,G0,HH,C,DX,SURNORMU
113 !
114  INTEGER I1,I2,I3,I4,I5,I6,IELEM
115 !
116  INTRINSIC sqrt
117 !
118 !**********************************************************************
119 !
120  IF(su%ELM.NE.41) THEN
121  WRITE(lu,1001) su%ELM
122 1001 FORMAT(1x,'MT04PP (BIEF) : TYPE OF U NOT IMPLEMENTED: ',i6)
123  CALL plante(1)
124  stop
125  ENDIF
126  IF(sv%ELM.NE.41) THEN
127  WRITE(lu,2001) sv%ELM
128 2001 FORMAT(1x,'MT04PP (BIEF) : TYPE OF V NOT IMPLEMENTED: ',i6)
129  CALL plante(1)
130  stop
131  ENDIF
132 ! HERE WSCONV WHICH IS IN FACT DEFINED PER LAYER
133  IF(sw%ELM.NE.41) THEN
134  WRITE(lu,3001) sw%ELM
135 3001 FORMAT(1x,'MT04PP (BIEF) : TYPE OF W NOT IMPLEMENTED: ',i6)
136  CALL plante(1)
137  stop
138  ENDIF
139 !
140  IF(formul(1:7).EQ.'MAUGUG2') THEN
141 !
142 ! LOOP ON THE 3D ELEMENTS
143 !
144  DO ielem = 1 , nelem
145 !
146  i1 = ikle(ielem,1)
147  i2 = ikle(ielem,2)
148  i3 = ikle(ielem,3)
149  i4 = ikle(ielem,4)
150  i5 = ikle(ielem,5)
151  i6 = ikle(ielem,6)
152 !
153 ! X2 = X(I2) - X(I1)
154 ! X3 = X(I3) - X(I1)
155 ! Y2 = Y(I2) - Y(I1)
156 ! Y3 = Y(I3) - Y(I1)
157  x2 = x(ielem,2)
158  x3 = x(ielem,3)
159  y2 = y(ielem,2)
160  y3 = y(ielem,3)
161 !
162  u0 = (u(i1)+u(i2)+u(i3)+u(i4)+u(i5)+u(i6))/6.d0
163  v0 = (v(i1)+v(i2)+v(i3)+v(i4)+v(i5)+v(i6))/6.d0
164 ! AVERAGE HEIGHT OF PRISM
165  hh = (z(i4)-z(i1)+z(i5)-z(i2)+z(i6)-z(i3))/3.d0
166  c = xmul*hh/24.d0/surfac(ielem)
167 !
168  t(ielem,1)=2*c*(u0*y3-u0*y2-v0*x3+v0*x2)**2
169  t(ielem,2)=2*c*(-u0*y3+v0*x3)**2
170  t(ielem,3)=2*c*(-u0*y2+v0*x2)**2
171  t(ielem,4)=t(ielem,1)
172  t(ielem,5)=t(ielem,2)
173  t(ielem,6)=t(ielem,3)
174 !
175  xm(ielem,01)= 2*(-u0*y3+v0*x3)*(u0*y3-u0*y2-v0*x3+v0*x2)*c
176  xm(ielem,02)=-2*(-u0*y2+v0*x2)*(u0*y3-u0*y2-v0*x3+v0*x2)*c
177  xm(ielem,03)=(u0*y3-u0*y2-v0*x3+v0*x2)*(u0*y3-u0*y2-v0*x3+v0*x2)*c
178  xm(ielem,04)=(-u0*y3+v0*x3)*(u0*y3-u0*y2-v0*x3+v0*x2)*c
179  xm(ielem,05)=-(-u0*y2+v0*x2)*(u0*y3-u0*y2-v0*x3+v0*x2)*c
180  xm(ielem,06)=-2*(-u0*y2+v0*x2)*(-u0*y3+v0*x3)*c
181  xm(ielem,07)=(u0*y3-u0*y2-v0*x3+v0*x2)*(-u0*y3+v0*x3)*c
182  xm(ielem,08)=(-u0*y3+v0*x3)*(-u0*y3+v0*x3)*c
183  xm(ielem,09)=-(-u0*y2+v0*x2)*(-u0*y3+v0*x3)*c
184  xm(ielem,10)=-(u0*y3-u0*y2-v0*x3+v0*x2)*(-u0*y2+v0*x2)*c
185  xm(ielem,11)=-(-u0*y3+v0*x3)*(-u0*y2+v0*x2)*c
186  xm(ielem,12)=(-u0*y2+v0*x2)*(-u0*y2+v0*x2)*c
187  xm(ielem,13)= 2*(-u0*y3+v0*x3)*(u0*y3-u0*y2-v0*x3+v0*x2)*c
188  xm(ielem,14)=-2*(-u0*y2+v0*x2)*(u0*y3-u0*y2-v0*x3+v0*x2)*c
189  xm(ielem,15)=-2*(-u0*y2+v0*x2)*(-u0*y3+v0*x3)*c
190 !
191  ENDDO
192 !
193  ELSEIF(formul(1:7).EQ.'MAUGUG1') THEN
194 !
195  DO ielem = 1 , nelem
196 !
197  i1 = ikle(ielem,1)
198  i2 = ikle(ielem,2)
199  i3 = ikle(ielem,3)
200  i4 = ikle(ielem,4)
201  i5 = ikle(ielem,5)
202  i6 = ikle(ielem,6)
203 !
204 ! X2 = X(I2) - X(I1)
205 ! X3 = X(I3) - X(I1)
206 ! Y2 = Y(I2) - Y(I1)
207 ! Y3 = Y(I3) - Y(I1)
208  x2 = x(ielem,2)
209  x3 = x(ielem,3)
210  y2 = y(ielem,2)
211  y3 = y(ielem,3)
212 !
213 ! VELOCITIES CONSIDERED CONSTANT
214  u0 = (u(i1)+u(i2)+u(i3)+u(i4)+u(i5)+u(i6))/6.d0
215  v0 = (v(i1)+v(i2)+v(i3)+v(i4)+v(i5)+v(i6))/6.d0
216 !
217  surnormu=1.d0/max(sqrt(u0**2+v0**2),1.d-8)
218  dx=sqrt(2.d0*surfac(ielem))
219 ! HERE F0 AND G0 MUST BE
220 ! DX*U/(2*NORME(U)) AND DY*V/(2*NORME(U)) BUT CONSIDERS DY=DX
221 ! APPROXIMATED BY SQUARE ROOT OF (2*TRIANGLE AREA)
222  f0 = 0.5d0*dx*u0*surnormu
223  g0 = 0.5d0*dx*v0*surnormu
224 ! AVERAGE HEIGHT OF PRISM
225  hh = (z(i4)-z(i1)+z(i5)-z(i2)+z(i6)-z(i3))/3.d0
226  c = xmul*hh/24.d0/surfac(ielem)
227 !
228  t(ielem,1)=2*(u0*y3-u0*y2-v0*x3+v0*x2)*(f0*y3-f0*y2-g0*x3+g0*x2)*c
229  t(ielem,2)=2*(-u0*y3+v0*x3)*(-f0*y3+g0*x3)*c
230  t(ielem,3)=2*(-u0*y2+v0*x2)*(-f0*y2+g0*x2)*c
231  t(ielem,4)=t(ielem,1)
232  t(ielem,5)=t(ielem,2)
233  t(ielem,6)=t(ielem,3)
234 !
235  xm(ielem,01)=2*(-u0*y3+v0*x3)*(f0*y3-f0*y2-g0*x3+g0*x2)*c
236  xm(ielem,02)=-2*(-u0*y2+v0*x2)*(f0*y3-f0*y2-g0*x3+g0*x2)*c
237  xm(ielem,03)=(u0*y3-u0*y2-v0*x3+v0*x2)*(f0*y3-f0*y2-g0*x3+g0*x2)*c
238  xm(ielem,04)=(-u0*y3+v0*x3)*(f0*y3-f0*y2-g0*x3+g0*x2)*c
239  xm(ielem,05)=-(-u0*y2+v0*x2)*(f0*y3-f0*y2-g0*x3+g0*x2)*c
240  xm(ielem,06)=-2*(-u0*y2+v0*x2)*(-f0*y3+g0*x3)*c
241  xm(ielem,07)=(u0*y3-u0*y2-v0*x3+v0*x2)*(-f0*y3+g0*x3)*c
242  xm(ielem,08)=(-u0*y3+v0*x3)*(-f0*y3+g0*x3)*c
243  xm(ielem,09)=-(-u0*y2+v0*x2)*(-f0*y3+g0*x3)*c
244  xm(ielem,10)=-(u0*y3-u0*y2-v0*x3+v0*x2)*(-f0*y2+g0*x2)*c
245  xm(ielem,11)=-(-u0*y3+v0*x3)*(-f0*y2+g0*x2)*c
246  xm(ielem,12)=(-u0*y2+v0*x2)*(-f0*y2+g0*x2)*c
247  xm(ielem,13)=2*(-u0*y3+v0*x3)*(f0*y3-f0*y2-g0*x3+g0*x2)*c
248  xm(ielem,14)=-2*(-u0*y2+v0*x2)*(f0*y3-f0*y2-g0*x3+g0*x2)*c
249  xm(ielem,15)=-2*(-u0*y2+v0*x2)*(-f0*y3+g0*x3)*c
250 !
251  xm(ielem,16)= 2*(u0*y3-u0*y2-v0*x3+v0*x2)*(-f0*y3+g0*x3)*c
252  xm(ielem,17)= -2*(u0*y3-u0*y2-v0*x3+v0*x2)*(-f0*y2+g0*x2)*c
253  xm(ielem,21)= -2*(-u0*y3+v0*x3)*(-f0*y2+g0*x2)*c
254  xm(ielem,18)=(u0*y3-u0*y2-v0*x3+v0*x2)*(f0*y3-f0*y2-g0*x3+g0*x2)*c
255  xm(ielem,22)= (-u0*y3+v0*x3)*(f0*y3-f0*y2-g0*x3+g0*x2)*c
256  xm(ielem,25)= -(-u0*y2+v0*x2)*(f0*y3-f0*y2-g0*x3+g0*x2)*c
257  xm(ielem,19)= (u0*y3-u0*y2-v0*x3+v0*x2)*(-f0*y3+g0*x3)*c
258  xm(ielem,23)= (-u0*y3+v0*x3)*(-f0*y3+g0*x3)*c
259  xm(ielem,26)= -(-u0*y2+v0*x2)*(-f0*y3+g0*x3)*c
260  xm(ielem,28)= 2*(u0*y3-u0*y2-v0*x3+v0*x2)*(-f0*y3+g0*x3)*c
261  xm(ielem,20)= -(u0*y3-u0*y2-v0*x3+v0*x2)*(-f0*y2+g0*x2)*c
262  xm(ielem,24)= -(-u0*y3+v0*x3)*(-f0*y2+g0*x2)*c
263  xm(ielem,27)= (-u0*y2+v0*x2)*(-f0*y2+g0*x2)*c
264  xm(ielem,29)= -2*(u0*y3-u0*y2-v0*x3+v0*x2)*(-f0*y2+g0*x2)*c
265  xm(ielem,30)= -2*(-u0*y3+v0*x3)*(-f0*y2+g0*x2)*c
266 !
267  ENDDO
268 !
269  ELSE
270  WRITE(lu,4001) formul
271 4001 FORMAT(1x,'MT04PP (BIEF) : UNEXPECTED FORMULA: ',a16)
272  CALL plante(1)
273  stop
274  ENDIF
275 !
276 !-----------------------------------------------------------------------
277 !
278  RETURN
279  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mt04pp
subroutine mt04pp(T, XM, XMUL, SU, SV, SW, U, V, X, Y, Z, SURFAC, IKLE, NELEM, NELMAX, FORMUL)
Definition: mt04pp.f:7
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:06 for The TELEMAC-MASCARET system by doxygen 1.8.13