The TELEMAC-MASCARET system  trunk
char_weak.F
Go to the documentation of this file.
1 ! ********************
2  SUBROUTINE char_weak
3 ! ********************
4 !
5  &(ftild,ftild_weak,surfac,ikle,npoin,nelem,nelmax,ng,ngauss,
6  & mesh,t2,tb,agglo,ielm,nplan,z,cv1,am1,slv,unsv,listin,solv)
7 !
8 !***********************************************************************
9 ! BIEF V6P3 21/08/2010
10 !***********************************************************************
11 !
12 !brief Completing the weak form of the method of characteristics
13 !+ after advection of the Gauss points.
14 !
15 !history J-M HERVOUET (EDF R&D, LNHE)
16 !+ 23/04/2013
17 !+ V6P3
18 !+ First version.
19 !
20 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
21 !| AM1 |-->| BIEF_OBJ WORK MATRIX
22 !| CV1 |-->| BIEF_OBJ WORK ARRAY
23 !| FTILD |-->| BLOCK OF RESULTS
24 !| FTILD_WEAK |-->| BLOCK OF RESULTS FOR ADVECTED GAUSS POINTS
25 !| IELM |-->| TYPE OF ELEMENT
26 !| IKLE |-->| CONNECTIVITY TABLE FOR ALL POINTS
27 !| NELEM |-->| NOMBRE D'ELEMENTS
28 !| NELMAX |-->| NOMBRE MAXIMUM D'ELEMENTS
29 !| NG |-->| TOTAL NUMBER OF GAUSS POINTS IN THE MESH
30 !| NGAUSS |-->| NUMBER OF GAUSS POINTS PER ELEMENT
31 !| NPOIN |-->| NUMBER OF POINTS IN THE MESH
32 !| SOLV |-->| IF YES, SOLVE THE LINEAR SYSTEM
33 !| SURFAC |-->| AREA OF ELEMENTS
34 !| Z |-->| ELEVATIONS OF POINTS IN THE MESH.
35 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
36 !
37  USE bief
38 !
39  USE declarations_special
40  IMPLICIT NONE
41 !
42 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
43 !
44  INTEGER, INTENT(IN) :: NELEM,NELMAX,NPOIN,NG,IELM
45  INTEGER, INTENT(IN) :: NPLAN,NGAUSS
46 ! HERE IKLE2 AND NELMAX2
47  INTEGER, INTENT(IN) :: IKLE(nelmax,3)
48  LOGICAL, INTENT(IN) :: LISTIN,SOLV
49  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax),AGGLO,Z(*)
50  TYPE(bief_obj), INTENT(IN) :: FTILD_WEAK,UNSV
51  TYPE(bief_obj), INTENT(INOUT) :: FTILD,T2,TB,CV1,AM1
52  TYPE(bief_mesh), INTENT(INOUT) :: MESH
53  TYPE(slvcfg), INTENT(INOUT) :: SLV
54 !
55 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
56 !
57  INTEGER IELEM,I1,I2,I3,I4,I5,I6,I,IG,IPLAN
58  DOUBLE PRECISION TIERS,A,B,C,D,H1,H2,H3,A1,A2,A3
59  DOUBLE PRECISION WEIGH1,WEIGH2,WEIGH3
60 !
61  tiers=1.d0/3.d0
62 !
63 !-----------------------------------------------------------------------
64 !
65 ! INITIALISATION
66 !
67  CALL cpstvc(ftild,cv1)
68 !
69  DO i=1,npoin*nplan
70  cv1%R(i)=0.d0
71  ENDDO
72 !
73  IF(ng.NE.nelem*ngauss.AND.ielm.EQ.11) THEN
74  WRITE(lu,*) 'CHAR_WEAK: BAD NUMBER OF POINTS'
75  WRITE(lu,*) 'NG=',ng,' NELEM=',nelem,' NGAUSS=',ngauss
76  CALL plante(1)
77  stop
78  ELSEIF(ng.NE.nelem*(nplan-1)*ngauss.AND.ielm.EQ.41) THEN
79  WRITE(lu,*) 'CHAR_WEAK: BAD NUMBER OF POINTS'
80  WRITE(lu,*) 'NG=',ng,' NELEM=',nelem,' NGAUSS=',ngauss
81  WRITE(lu,*) 'NPLAN=',nplan
82  CALL plante(1)
83  stop
84  ENDIF
85 !
86  IF(ngauss.EQ.1.AND.ielm.EQ.11) THEN
87 !
88 ! ASSEMBLING (3 BASES PER ELEMENT)
89 ! HERE THE VALUE OF THE BASIS IS 1/3 FOR ALL 3 OF THEM
90 ! AND THE WEIGHTS ARE ALL SURFAC
91 !
92  DO ielem=1,nelem
93  i1=ikle(ielem,1)
94  i2=ikle(ielem,2)
95  i3=ikle(ielem,3)
96  cv1%R(i1)=cv1%R(i1)+surfac(ielem)*tiers*ftild_weak%R(ielem)
97  cv1%R(i2)=cv1%R(i2)+surfac(ielem)*tiers*ftild_weak%R(ielem)
98  cv1%R(i3)=cv1%R(i3)+surfac(ielem)*tiers*ftild_weak%R(ielem)
99  ENDDO
100 !
101  ELSEIF(ngauss.EQ.3.AND.ielm.EQ.11) THEN
102 !
103 ! ASSEMBLING (3 BASES PER ELEMENT)
104 ! HERE THE WEIGHTS ARE ALL SURFAC/3
105 ! THE VALUES IF THE TEST FUNCTIONS AT GAUSS POINTS ARE ON THE RIGHT
106 !
107  ig=0
108  DO ielem=1,nelem
109  i1=ikle(ielem,1)
110  i2=ikle(ielem,2)
111  i3=ikle(ielem,3)
112 ! POINT DE GAUSS 1
113  ig=ig+1
114  cv1%R(i1)=
115  & cv1%R(i1)+surfac(ielem)*tiers*ftild_weak%R(ig)*2.d0/3.d0
116  cv1%R(i2)=
117  & cv1%R(i2)+surfac(ielem)*tiers*ftild_weak%R(ig)*1.d0/6.d0
118  cv1%R(i3)=
119  & cv1%R(i3)+surfac(ielem)*tiers*ftild_weak%R(ig)*1.d0/6.d0
120 ! POINT DE GAUSS 2
121  ig=ig+1
122  cv1%R(i1)=
123  & cv1%R(i1)+surfac(ielem)*tiers*ftild_weak%R(ig)*1.d0/6.d0
124  cv1%R(i2)=
125  & cv1%R(i2)+surfac(ielem)*tiers*ftild_weak%R(ig)*2.d0/3.d0
126  cv1%R(i3)=
127  & cv1%R(i3)+surfac(ielem)*tiers*ftild_weak%R(ig)*1.d0/6.d0
128 ! POINT DE GAUSS 3
129  ig=ig+1
130  cv1%R(i1)=
131  & cv1%R(i1)+surfac(ielem)*tiers*ftild_weak%R(ig)*1.d0/6.d0
132  cv1%R(i2)=
133  & cv1%R(i2)+surfac(ielem)*tiers*ftild_weak%R(ig)*1.d0/6.d0
134  cv1%R(i3)=
135  & cv1%R(i3)+surfac(ielem)*tiers*ftild_weak%R(ig)*2.d0/3.d0
136  ENDDO
137 !
138  ELSEIF(ngauss.EQ.4.AND.ielm.EQ.11) THEN
139 !
140 ! ASSEMBLING (3 BASES PER ELEMENT)
141 ! THE WEIGHS ARE:
142 ! -27*SURFAC/48 FOR POINT 1
143 ! 25*SURFAC/48 FOR POINT 2,3 AND 4
144 !
145  ig=0
146  DO ielem=1,nelem
147  i1=ikle(ielem,1)
148  i2=ikle(ielem,2)
149  i3=ikle(ielem,3)
150 ! MONOTONICITY IF FORMULA APPROXIMATE ?
151  weigh1=-27.d0*surfac(ielem)/48.d0
152  weigh2= 25.d0*surfac(ielem)/48.d0
153 ! POINT DE GAUSS 1
154  ig=ig+1
155  cv1%R(i1)=cv1%R(i1)+weigh1*ftild_weak%R(ig)*1.d0/3.d0
156  cv1%R(i2)=cv1%R(i2)+weigh1*ftild_weak%R(ig)*1.d0/3.d0
157  cv1%R(i3)=cv1%R(i3)+weigh1*ftild_weak%R(ig)*1.d0/3.d0
158 ! POINT DE GAUSS 2
159  ig=ig+1
160  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*3.d0/5.d0
161  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*1.d0/5.d0
162  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*1.d0/5.d0
163 ! POINT DE GAUSS 3
164  ig=ig+1
165  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*1.d0/5.d0
166  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*3.d0/5.d0
167  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*1.d0/5.d0
168 ! POINT DE GAUSS 4
169  ig=ig+1
170  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*1.d0/5.d0
171  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*1.d0/5.d0
172  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*3.d0/5.d0
173  ENDDO
174 !
175  ELSEIF(ngauss.EQ.6.AND.ielm.EQ.11) THEN
176 !
177 ! ASSEMBLING (3 BASES PER ELEMENT)
178 ! THE WEIGHS ARE:
179 ! WEIGH1 FOR POINT 1,2,3
180 ! WEIGH2 FOR POINT 4,5,6
181 !
182  ig=0
183  a=0.445948490915965d0
184  b=0.091576213509771d0
185  DO ielem=1,nelem
186  i1=ikle(ielem,1)
187  i2=ikle(ielem,2)
188  i3=ikle(ielem,3)
189  weigh1=surfac(ielem)*2.d0*0.111690794839005d0
190  weigh2=surfac(ielem)*2.d0*0.054975871827661d0
191 ! POINT DE GAUSS 1
192  ig=ig+1
193  cv1%R(i1)=cv1%R(i1)+weigh1*ftild_weak%R(ig)*(1.d0-a-a)
194  cv1%R(i2)=cv1%R(i2)+weigh1*ftild_weak%R(ig)*a
195  cv1%R(i3)=cv1%R(i3)+weigh1*ftild_weak%R(ig)*a
196 ! POINT DE GAUSS 2
197  ig=ig+1
198  cv1%R(i1)=cv1%R(i1)+weigh1*ftild_weak%R(ig)*a
199  cv1%R(i2)=cv1%R(i2)+weigh1*ftild_weak%R(ig)*(1.d0-a-a)
200  cv1%R(i3)=cv1%R(i3)+weigh1*ftild_weak%R(ig)*a
201 ! POINT DE GAUSS 3
202  ig=ig+1
203  cv1%R(i1)=cv1%R(i1)+weigh1*ftild_weak%R(ig)*a
204  cv1%R(i2)=cv1%R(i2)+weigh1*ftild_weak%R(ig)*a
205  cv1%R(i3)=cv1%R(i3)+weigh1*ftild_weak%R(ig)*(1.d0-a-a)
206 ! POINT DE GAUSS 4
207  ig=ig+1
208  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*(1.d0-b-b)
209  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*b
210  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*b
211 ! POINT DE GAUSS 5
212  ig=ig+1
213  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*b
214  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*(1.d0-b-b)
215  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*b
216 ! POINT DE GAUSS 6
217  ig=ig+1
218  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*b
219  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*b
220  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*(1.d0-b-b)
221  ENDDO
222 !
223  ELSEIF(ngauss.EQ.7.AND.ielm.EQ.11) THEN
224 !
225 ! ASSEMBLING (3 BASES PER ELEMENT)
226 ! THE WEIGHS ARE:
227 ! WEIGH1 FOR POINT 1,2,3
228 ! WEIGH2 FOR POINT 4,5,6
229 !
230  ig=0
231  a=(6.d0+sqrt(15.d0))/21.d0
232  b=4.d0/7.d0-a
233  a1=9.d0/80.d0
234  a2=(155.d0+sqrt(15.d0))/2400.d0
235  a3=31.d0/240.d0-a2
236  DO ielem=1,nelem
237  i1=ikle(ielem,1)
238  i2=ikle(ielem,2)
239  i3=ikle(ielem,3)
240  weigh1=surfac(ielem)*2.d0*a1
241  weigh2=surfac(ielem)*2.d0*a2
242  weigh3=surfac(ielem)*2.d0*a3
243 ! POINT DE GAUSS 1
244  ig=ig+1
245  cv1%R(i1)=cv1%R(i1)+weigh1*ftild_weak%R(ig)*tiers
246  cv1%R(i2)=cv1%R(i2)+weigh1*ftild_weak%R(ig)*tiers
247  cv1%R(i3)=cv1%R(i3)+weigh1*ftild_weak%R(ig)*tiers
248 ! POINT DE GAUSS 2
249  ig=ig+1
250  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*(1.d0-a-a)
251  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*a
252  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*a
253 ! POINT DE GAUSS 3
254  ig=ig+1
255  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*a
256  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*(1.d0-a-a)
257  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*a
258 ! POINT DE GAUSS 4
259  ig=ig+1
260  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*a
261  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*a
262  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*(1.d0-a-a)
263 ! POINT DE GAUSS 5
264  ig=ig+1
265  cv1%R(i1)=cv1%R(i1)+weigh3*ftild_weak%R(ig)*(1.d0-b-b)
266  cv1%R(i2)=cv1%R(i2)+weigh3*ftild_weak%R(ig)*b
267  cv1%R(i3)=cv1%R(i3)+weigh3*ftild_weak%R(ig)*b
268 ! POINT DE GAUSS 6
269  ig=ig+1
270  cv1%R(i1)=cv1%R(i1)+weigh3*ftild_weak%R(ig)*b
271  cv1%R(i2)=cv1%R(i2)+weigh3*ftild_weak%R(ig)*(1.d0-b-b)
272  cv1%R(i3)=cv1%R(i3)+weigh3*ftild_weak%R(ig)*b
273 ! POINT DE GAUSS 7
274  ig=ig+1
275  cv1%R(i1)=cv1%R(i1)+weigh3*ftild_weak%R(ig)*b
276  cv1%R(i2)=cv1%R(i2)+weigh3*ftild_weak%R(ig)*b
277  cv1%R(i3)=cv1%R(i3)+weigh3*ftild_weak%R(ig)*(1.d0-b-b)
278  ENDDO
279 !
280  ELSEIF(ngauss.EQ.12.AND.ielm.EQ.11) THEN
281 !
282  ig=0
283  a=0.063089014491502d0
284  b=0.249286745170910d0
285  c=0.310352451033785d0
286  d=0.053145049844816d0
287  a1=0.025422453185103d0
288  a2=0.058393137863189d0
289  a3=0.041425537809187d0
290  a3=(1.d0-6.d0*a1-6.d0*a2)/12.d0
291  DO ielem=1,nelem
292  i1=ikle(ielem,1)
293  i2=ikle(ielem,2)
294  i3=ikle(ielem,3)
295  weigh1=surfac(ielem)*2.d0*a1
296  weigh2=surfac(ielem)*2.d0*a2
297  weigh3=surfac(ielem)*2.d0*a3
298 ! POINT DE GAUSS 1
299  ig=ig+1
300  cv1%R(i1)=cv1%R(i1)+weigh1*ftild_weak%R(ig)*(1.d0-a-a)
301  cv1%R(i2)=cv1%R(i2)+weigh1*ftild_weak%R(ig)*a
302  cv1%R(i3)=cv1%R(i3)+weigh1*ftild_weak%R(ig)*a
303 ! POINT DE GAUSS 2
304  ig=ig+1
305  cv1%R(i1)=cv1%R(i1)+weigh1*ftild_weak%R(ig)*a
306  cv1%R(i2)=cv1%R(i2)+weigh1*ftild_weak%R(ig)*(1.d0-a-a)
307  cv1%R(i3)=cv1%R(i3)+weigh1*ftild_weak%R(ig)*a
308 ! POINT DE GAUSS 3
309  ig=ig+1
310  cv1%R(i1)=cv1%R(i1)+weigh1*ftild_weak%R(ig)*a
311  cv1%R(i2)=cv1%R(i2)+weigh1*ftild_weak%R(ig)*a
312  cv1%R(i3)=cv1%R(i3)+weigh1*ftild_weak%R(ig)*(1.d0-a-a)
313 ! POINT DE GAUSS 4
314  ig=ig+1
315  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*(1.d0-b-b)
316  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*b
317  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*b
318 ! POINT DE GAUSS 5
319  ig=ig+1
320  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*b
321  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*(1.d0-b-b)
322  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*b
323 ! POINT DE GAUSS 6
324  ig=ig+1
325  cv1%R(i1)=cv1%R(i1)+weigh2*ftild_weak%R(ig)*b
326  cv1%R(i2)=cv1%R(i2)+weigh2*ftild_weak%R(ig)*b
327  cv1%R(i3)=cv1%R(i3)+weigh2*ftild_weak%R(ig)*(1.d0-b-b)
328 ! POINT DE GAUSS 7
329  ig=ig+1
330  cv1%R(i1)=cv1%R(i1)+weigh3*ftild_weak%R(ig)*(1.d0-c-d)
331  cv1%R(i2)=cv1%R(i2)+weigh3*ftild_weak%R(ig)*c
332  cv1%R(i3)=cv1%R(i3)+weigh3*ftild_weak%R(ig)*d
333 ! POINT DE GAUSS 8
334  ig=ig+1
335  cv1%R(i1)=cv1%R(i1)+weigh3*ftild_weak%R(ig)*(1.d0-c-d)
336  cv1%R(i2)=cv1%R(i2)+weigh3*ftild_weak%R(ig)*d
337  cv1%R(i3)=cv1%R(i3)+weigh3*ftild_weak%R(ig)*c
338 ! POINT DE GAUSS 9
339  ig=ig+1
340  cv1%R(i1)=cv1%R(i1)+weigh3*ftild_weak%R(ig)*d
341  cv1%R(i2)=cv1%R(i2)+weigh3*ftild_weak%R(ig)*(1.d0-c-d)
342  cv1%R(i3)=cv1%R(i3)+weigh3*ftild_weak%R(ig)*c
343 ! POINT DE GAUSS 10
344  ig=ig+1
345  cv1%R(i1)=cv1%R(i1)+weigh3*ftild_weak%R(ig)*c
346  cv1%R(i2)=cv1%R(i2)+weigh3*ftild_weak%R(ig)*(1.d0-c-d)
347  cv1%R(i3)=cv1%R(i3)+weigh3*ftild_weak%R(ig)*d
348 ! POINT DE GAUSS 11
349  ig=ig+1
350  cv1%R(i1)=cv1%R(i1)+weigh3*ftild_weak%R(ig)*d
351  cv1%R(i2)=cv1%R(i2)+weigh3*ftild_weak%R(ig)*c
352  cv1%R(i3)=cv1%R(i3)+weigh3*ftild_weak%R(ig)*(1.d0-c-d)
353 ! POINT DE GAUSS 12
354  ig=ig+1
355  cv1%R(i1)=cv1%R(i1)+weigh3*ftild_weak%R(ig)*c
356  cv1%R(i2)=cv1%R(i2)+weigh3*ftild_weak%R(ig)*d
357  cv1%R(i3)=cv1%R(i3)+weigh3*ftild_weak%R(ig)*(1.d0-c-d)
358  ENDDO
359 !
360  ELSEIF(ngauss.EQ.6.AND.ielm.EQ.41) THEN
361 !
362 ! ASSEMBLING (3 BASES PER ELEMENT)
363 ! THE WEIGHS ARE:
364 ! WEIGH1 FOR POINT 1,2,3
365 ! WEIGH2 FOR POINT 4,5,6
366 !
367  a=1.d0/6.d0
368  b=2.d0/3.d0
369  c=0.5d0*(1.d0-1.d0/sqrt(3.d0))
370  d=0.5d0*(1.d0+1.d0/sqrt(3.d0))
371 !
372  ig=0
373  DO iplan=1,nplan-1
374  DO ielem=1,nelem
375  i1=ikle(ielem,1)+(iplan-1)*npoin
376  i2=ikle(ielem,2)+(iplan-1)*npoin
377  i3=ikle(ielem,3)+(iplan-1)*npoin
378  i4=i1+npoin
379  i5=i2+npoin
380  i6=i3+npoin
381  h1=z(i4)-z(i1)
382  h2=z(i5)-z(i2)
383  h3=z(i6)-z(i3)
384  weigh1=surfac(ielem)*(h1*(1.d0-a-a)+h2*a+h3*a)/6.d0
385  weigh2=surfac(ielem)*(h1*(1.d0-b-a)+h2*b+h3*a)/6.d0
386  weigh3=surfac(ielem)*(h1*(1.d0-a-b)+h2*a+h3*b)/6.d0
387 ! POINT DE GAUSS 1
388  ig=ig+1
389  cv1%R(i1)=cv1%R(i1)
390  & +weigh1*ftild_weak%R(ig)*(1.d0-c)*(1.d0-a-a)
391  cv1%R(i2)=cv1%R(i2)
392  & +weigh1*ftild_weak%R(ig)*(1.d0-c)* a
393  cv1%R(i3)=cv1%R(i3)
394  & +weigh1*ftild_weak%R(ig)*(1.d0-c)* a
395  cv1%R(i4)=cv1%R(i4)
396  & +weigh1*ftild_weak%R(ig)* c *(1.d0-a-a)
397  cv1%R(i5)=cv1%R(i5)
398  & +weigh1*ftild_weak%R(ig)* c * a
399  cv1%R(i6)=cv1%R(i6)
400  & +weigh1*ftild_weak%R(ig)* c * a
401 ! POINT DE GAUSS 2
402  ig=ig+1
403  cv1%R(i1)=cv1%R(i1)
404  & +weigh2*ftild_weak%R(ig)*(1.d0-c)*(1.d0-b-a)
405  cv1%R(i2)=cv1%R(i2)
406  & +weigh2*ftild_weak%R(ig)*(1.d0-c)* b
407  cv1%R(i3)=cv1%R(i3)
408  & +weigh2*ftild_weak%R(ig)*(1.d0-c)* a
409  cv1%R(i4)=cv1%R(i4)
410  & +weigh2*ftild_weak%R(ig)* c *(1.d0-b-a)
411  cv1%R(i5)=cv1%R(i5)
412  & +weigh2*ftild_weak%R(ig)* c * b
413  cv1%R(i6)=cv1%R(i6)
414  & +weigh2*ftild_weak%R(ig)* c * a
415 ! POINT DE GAUSS 3
416  ig=ig+1
417  cv1%R(i1)=cv1%R(i1)
418  & +weigh3*ftild_weak%R(ig)*(1.d0-c)*(1.d0-a-b)
419  cv1%R(i2)=cv1%R(i2)
420  & +weigh3*ftild_weak%R(ig)*(1.d0-c)* a
421  cv1%R(i3)=cv1%R(i3)
422  & +weigh3*ftild_weak%R(ig)*(1.d0-c)* b
423  cv1%R(i4)=cv1%R(i4)
424  & +weigh3*ftild_weak%R(ig)* c *(1.d0-a-b)
425  cv1%R(i5)=cv1%R(i5)
426  & +weigh3*ftild_weak%R(ig)* c * a
427  cv1%R(i6)=cv1%R(i6)
428  & +weigh3*ftild_weak%R(ig)* c * b
429 ! POINT DE GAUSS 4
430  ig=ig+1
431  cv1%R(i1)=cv1%R(i1)
432  & +weigh1*ftild_weak%R(ig)*(1.d0-d)*(1.d0-a-a)
433  cv1%R(i2)=cv1%R(i2)
434  & +weigh1*ftild_weak%R(ig)*(1.d0-d)* a
435  cv1%R(i3)=cv1%R(i3)
436  & +weigh1*ftild_weak%R(ig)*(1.d0-d)* a
437  cv1%R(i4)=cv1%R(i4)
438  & +weigh1*ftild_weak%R(ig)* d *(1.d0-a-a)
439  cv1%R(i5)=cv1%R(i5)
440  & +weigh1*ftild_weak%R(ig)* d * a
441  cv1%R(i6)=cv1%R(i6)
442  & +weigh1*ftild_weak%R(ig)* d * a
443 ! POINT DE GAUSS 5
444  ig=ig+1
445  cv1%R(i1)=cv1%R(i1)
446  & +weigh2*ftild_weak%R(ig)*(1.d0-d)*(1.d0-b-a)
447  cv1%R(i2)=cv1%R(i2)
448  & +weigh2*ftild_weak%R(ig)*(1.d0-d)* b
449  cv1%R(i3)=cv1%R(i3)
450  & +weigh2*ftild_weak%R(ig)*(1.d0-d)* a
451  cv1%R(i4)=cv1%R(i4)
452  & +weigh2*ftild_weak%R(ig)* d *(1.d0-b-a)
453  cv1%R(i5)=cv1%R(i5)
454  & +weigh2*ftild_weak%R(ig)* d * b
455  cv1%R(i6)=cv1%R(i6)
456  & +weigh2*ftild_weak%R(ig)* d * a
457 ! POINT DE GAUSS 6
458  ig=ig+1
459  cv1%R(i1)=cv1%R(i1)
460  & +weigh3*ftild_weak%R(ig)*(1.d0-d)*(1.d0-a-b)
461  cv1%R(i2)=cv1%R(i2)
462  & +weigh3*ftild_weak%R(ig)*(1.d0-d)* a
463  cv1%R(i3)=cv1%R(i3)
464  & +weigh3*ftild_weak%R(ig)*(1.d0-d)* b
465  cv1%R(i4)=cv1%R(i4)
466  & +weigh3*ftild_weak%R(ig)* d *(1.d0-a-b)
467  cv1%R(i5)=cv1%R(i5)
468  & +weigh3*ftild_weak%R(ig)* d * a
469  cv1%R(i6)=cv1%R(i6)
470  & +weigh3*ftild_weak%R(ig)* d * b
471  ENDDO
472  ENDDO
473 !
474  ELSE
475 !
476  WRITE(lu,11) ngauss
477 11 FORMAT(1x,'CHAR_WEAK: OPTION NOT IMPLEMENTED:',i6)
478  CALL plante(1)
479  stop
480 !
481  ENDIF
482 !
483 !-----------------------------------------------------------------------
484 !
485  IF(solv) THEN
486 !
487 ! MASS-MATRIX
488 !
489  CALL matrix(am1,'M=N ','MATMAS ',ielm,ielm,
490  & 1.d0,ftild,ftild,ftild,ftild,ftild,ftild,
491  & mesh,.false.,ftild)
492 !
493 ! PARTIALLY LUMPING THE MASS-MATRIX
494 !
495  CALL lump(t2,am1,mesh,agglo)
496  CALL om('M=CM ', m=am1, c=1.d0-agglo, mesh=mesh)
497  CALL om('M=M+D ', m=am1, d=t2, mesh=mesh)
498 !
499 ! INITIAL GUESS (AS IF MATRIX LUMPED...)
500 !
501  CALL os('X=Y ',x=ftild,y=cv1)
502  IF(ncsize.GT.1) CALL parcom(ftild,2,mesh)
503  CALL os('X=XY ',x=ftild,y=unsv)
504 !
505 ! SOLUTION OF THE SYSTEM (WITHOUT BOUNDARY CONDITIONS...)
506 !
507 #if defined COMPAD
508  CALL ad_solve(ftild,am1,cv1,tb,slv,listin,mesh,mesh%M)
509 #else
510  CALL solve(ftild,am1,cv1,tb,slv,listin,mesh,mesh%M)
511 #endif
512 !
513  ENDIF
514 !
515 !-----------------------------------------------------------------------
516 !
517  RETURN
518  END
ad_solve
subroutine ad_solve(X, A, B, TB, CFG, INFOGR, MESH, AUX)
Definition: ad_solve.F:33
declarations_special
Definition: declarations_special.F:3
solve
subroutine solve(X, A, B, TB, CFG, INFOGR, MESH, AUX)
Definition: solve.f:7
char_weak
subroutine char_weak(FTILD, FTILD_WEAK, SURFAC, IKLE, NPOIN, NELEM, NELMAX, NG, NGAUSS, MESH, T2, TB, AGGLO, IELM, NPLAN, Z, CV1, AM1, SLV, UNSV, LISTIN, SOLV)
Definition: char_weak.F:8
om
subroutine om(OP, M, N, D, C, MESH)
Definition: om.f:7
declarations_special::lu
integer lu
Definition: declarations_special.F:45
matrix
subroutine matrix(M, OP, FORMUL, IELM1, IELM2, XMUL, F, G, H, U, V, W, MESH, MSK, MASKEL)
Definition: matrix.f:7
cpstvc
subroutine cpstvc(X, Y)
Definition: cpstvc.f:7
os
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
parcom
subroutine parcom(X, ICOM, MESH)
Definition: parcom.f:7
lump
subroutine lump(DIAG, A, MESH, XMUL)
Definition: lump.f:7
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:36:31 for The TELEMAC-MASCARET system by doxygen 1.8.13