1 : C RLENCO SOURCE CHAT 05/01/13 03:01:24 5004
2 : SUBROUTINE RLENCO(MELSOM,MELCEN,MELLI1,MELLI2,INORM,MLEPOI,MLECOE)
3 : C************************************************************************
4 : C
5 : C PROJET : CASTEM 2000
6 : C
7 : C NOM : RLENCO
8 : C
9 : C DESCRIPTION : Appelle par GRADI2
10 : C
11 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
12 : C
13 : C AUTEUR : A. BECCANTINI
14 : C
15 : C************************************************************************
16 : C
17 : C This subroutine computes the coefficients which allow to perform
18 : C a linear exact reconstruction of each 'sommet' (vertex) value
19 : C Indeed, given the i-th sommet (V=MELSOM.NUM(1,i)), given the list
20 : C of its neighbors MLEPOI.LECT(i), we have to compute the matrix of
21 : C coefficients a(i,j) such that
22 : C
23 : C VAL(MELSOM.NUM(1,i)) = \sum_{neig1} a(i,neig1) * VAL(neig1) +
24 : C \sum_{neig2} a(i,neig2) * VALNOR(neig2)
25 : C
26 : C where neig1 is a 'CENTRE' point or a 'Dirichlet boundary condition'
27 : C point, neig2 is a 'Von Neumann boundary condition' point, VAL(j)
28 : C is the value of the function in the j-th neighbor, VALNOR(j) is
29 : C the value of grad(VAL).n is the j-th neighbor.
30 : C
31 : C
32 : C To compute these coefficients, we impose that VAL is linear with
33 : C respect to the coordinates of V and its neighbor. Then, since
34 : C there are less unknowns than equations, we use a least square
35 : C method. If the neighbor is a 'CENTRE' point or a 'Dirichlet
36 : C boundary condition' point, we impose that
37 : C
38 : C A + B (x_{neig} - x_V) + C (y_{neig} - y_V) + D (z_{neig} - z_V)
39 : C = VAL(neig)
40 : C
41 : C If the neighbor is a 'von Neumann boundary condition' point,
42 : C we impose that
43 : C
44 : C |r_{neig} - r_V| (B nx_{neig} + C ny_{neig} + D nz_{neig}) =
45 : C VALNOR(neig) * |r_{neig} - r_V|
46 : C
47 : c r = (x,y,z)
48 : c n = (nx, ny, nz) is the normal to the surface
49 : C
50 : C To determine A, B, C, D we have to solve
51 : C
52 : C MATA . tran(A, B, C, D) = tran(VAL(1), VAL(2), ... , VALNOR(j) /
53 : C (r_{j} - r_V)
54 : C
55 : C with MATA(1,*) = (1,x_1-x_V,y_1-y_V,z_1-z_V)
56 : C MATA(2,*) = (1,x_2-x_V,y_2-y_V,z_2-z_V)
57 : C ...
58 : C MATA(j,*) = (0,nx_{j}, ny_{j}, nz_{j})
59 : C ...
60 : C
61 : C To determine a(i,neig) we have to solve a linear system (for each
62 : C neig). If the neighbor is a 'CENTRE' point or a 'Dirichlet boundary
63 : C condition' point we solve
64 : C
65 : C MATA . tran(a(i,neig), b(i,neig), c(i,neig), d(i,neig)) = e_{neig}
66 : C
67 : C with e_{neig} = tran (0,0,...,0,1,0,...) (1 in the neig position)
68 : C
69 : C If the neighbor is a 'von Neumann boundary condition' points,
70 : C we solve
71 : C
72 : C MATA . tran(a(i,neig), b(i,neig), c(i,neig), d(i,neig)) =
73 : C e_{neig} (r_{neig} - r_V)
74 : C
75 : C Since each linear system involves MATA, we have to "pseudoinvert"
76 : C it.
77 : C
78 : C**** Inputs:
79 : C
80 : C MELSOM = the 'SOMMET' meleme
81 : C MELCEN = the 'CENTRE' meleme
82 : C MELLI1 = the support of Dirichlet boundary conditions
83 : C MELLI2 = the support of von Neumann boundary conditions
84 : C INOMR = CHAMPOINT des normales aux faces
85 : C MLEPOI = list of integers.
86 : C MLEPOI.LECT(i) points to the list neighbors of
87 : C MELSOM.NUM(1,i)
88 : C
89 : C**** Output:
90 : C
91 : C MLECOE = list of integers.
92 : C
93 : C MLECOE.LECT(i) points to the MLREEL MLRCOE which
94 : C contain the a(i,neig)
95 : C
96 : C**** Variables de COOPTIO
97 : C
98 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
99 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
100 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
101 : C & ,IECHO, IIMPI, IOSPI
102 : C & ,IDIM
103 : CC & ,MCOORD
104 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
105 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
106 : C & ,NORINC,NORVAL,NORIND,NORVAD
107 : C & ,NUCROU, IPSAUV, IFICLE, IPREFI
108 : C
109 : IMPLICIT INTEGER(I-N)
110 : INTEGER INORM, NSOMM, INOEU, NGS, IPCOOR, IVOI, NVOI, NGV, NLV
111 : & ,NLV1,NVMAX,NVMIN,IERSVD,IERR0,NLN
112 : & ,I1,I2,I3,I4
113 : REAL*8 XS,YS,ZS,XV,YV,ZV,DX,ERRTOL,SMAX,SMIN
114 : PARAMETER(ERRTOL=1.0D-16)
115 : CHARACTER*8 TYPE
116 : LOGICAL LOGSVD
117 : -INC CCOPTIO
118 : -INC SMCOORD
119 : -INC SMELEME
120 : -INC SMLREEL
121 : -INC SMLENTI
122 : -INC SMCHPOI
123 : INTEGER JG
124 : INTEGER N1,N2
125 : SEGMENT MATRIX
126 : REAL*8 MAT(N1,N2)
127 : ENDSEGMENT
128 : POINTEUR MELSOM.MELEME, MLEPOI.MLENTI,MELLI1.MELEME,MELLI2.MELEME
129 : & ,MLELI1.MLENTI,MLELI2.MLENTI,MELCEN.MELEME,MLECEN.MLENTI
130 : & ,MPNORM.MPOVAL, MATA.MATRIX,MATU.MATRIX,MATV.MATRIX
131 : & ,MLRCOE.MLREEL,MLRSIG.MLREEL,MLRTRA.MLREEL
132 : & ,MTAA.MATRIX, MINTAA.MATRIX,MLRIN1.MLREEL,MLRIN2.MLREEL
133 : & ,MLRIN3.MLREEL, MLECOE.MLENTI, MLRB.MLREEL
134 : & ,MELFAC.MELEME,MLEFAC.MLENTI
135 : C
136 : C**** We create the MLENTI for the centers
137 : C
137 : CALL KRIPAD(MELCEN,MLECEN)
138 : IF(IERR .NE. 0)GOTO 9999
139 : C En KRIPAD
140 : C SEGACT MELCEN
141 : C SEGINI MLECEN
142 : C
143 : C**** We create the MLENTI for the BC
144 : C
144 : CALL KRIPAD(MELLI1,MLELI1)
145 : IF(IERR .NE. 0)GOTO 9999
145 : CALL KRIPAD(MELLI2,MLELI2)
146 : IF(IERR .NE. 0)GOTO 9999
147 : C SEGACT MELLI1
148 : C SEGINI MLELI1
149 : C SEGACT MELLI2
150 : C SEGINI MLELI2
151 : C
152 : SEGACT MELSOM
153 : NSOMM=MELSOM.NUM(/2)
154 : JG=NSOMM
155 : SEGINI MLECOE
156 : SEGACT MLEPOI
157 : C
158 : C**** Le MPOVAL des normales
159 : C
159 : CALL LICHT(INORM,MPNORM,TYPE,MELFAC)
160 : C SEGACT MPNORM
161 : C
162 : C**** We create the MLENTI for the MELFAC
163 : C
163 : CALL KRIPAD(MELFAC,MLEFAC)
164 : IF(IERR .NE. 0)GOTO 9999
165 : C SEGACT MELFAC
166 : C SEGINI MLEFAC
167 : C
168 : C**** Loop on the sommet to compute NVMAX
169 : C (maximum number of neighbors)
170 : C
171 : NVMAX=0
172 : DO INOEU=1,NSOMM,1
173 : C
174 : C******* The neighbors of NGS
175 : C
176 : MLENTI=MLEPOI.LECT(INOEU)
177 : SEGACT MLENTI
178 : NVOI=MLENTI.LECT(/1)
179 : NVMAX=MAX(NVMAX,NVOI)
180 : ENDDO
181 : C
182 : C MATA = matrix to "pseudoinvert" (NVMAX,IDIM+1)
183 : C MATU = matrix of the singular right eigenvectors of MATA
184 : C (NVMAX,IDIM+1)
185 : C MATV = matrix of the singular left eigenvectors of MATA
186 : C (IDIM+1,IDIM+1)
187 : C But in invsvd.eso, MATV dimensions are NVMAX,IDIM+1
188 : C MLRSIG = singular values of MATA (IDIM+1)
189 : C MLRB = vector (NVMAX)
190 : C MLRB.PROG(j) = 1 if j is a center point or a 'Dirichlet
191 : C boundary condition' point
192 : C MLRB.PROG(j) = (r_j - r_V) id j is a 'von Neumann
193 : C boundary condition point.
194 : C
195 : C N.B. MATA = MATU MATSIG t(MATV)
196 : C If MATA is non singular,
197 : C inv(MATA) = MATV inv(MATSIG) t(MATU)
198 : C
199 : C MLRTRA temporary vector of invsvd.eso
200 : C NVMIN = IDIM + 1 (the most little dimension of the matrices)
201 : C
202 : N1=NVMAX
203 : N2=IDIM+1
204 : SEGINI MATA
205 : SEGINI MATU
206 : SEGINI MATV
207 : JG=NVMAX
208 : SEGINI MLRB
209 : JG=IDIM+1
210 : SEGINI MLRSIG
211 : SEGINI MLRTRA
212 : NVMIN=N2
213 : C
214 : C MTAA : [t(MATA).MATA]
215 : C MINTAA : [inve(t(MATA) MATA)]
216 : C MLRIN1,2,3 : "temporary vectors"
217 : C
218 : N1=NVMIN
219 : N2=NVMIN
220 : SEGINI MTAA
221 : SEGINI MINTAA
222 : JG=NVMIN
223 : SEGINI MLRIN1
224 : SEGINI MLRIN2
225 : SEGINI MLRIN3
226 : C
227 : C**** Loop on the sommet to compute the coefficient
228 : C
229 :
230 : DO INOEU=1,NSOMM,1
231 : NGS=MELSOM.NUM(1,INOEU)
232 : IPCOOR=((NGS-1)*(IDIM+1))+1
233 : XS=MCOORD.XCOOR(IPCOOR)
234 : YS=MCOORD.XCOOR(IPCOOR+1)
235 : IF(IDIM .EQ. 3) ZS=MCOORD.XCOOR(IPCOOR+2)
236 : C
237 : C******* The neighbors of NGS
238 : C
239 : MLENTI=MLEPOI.LECT(INOEU)
240 : NVOI=MLENTI.LECT(/1)
241 : C
242 : C******* The matrix where we store the coefficients
243 : C
244 : JG=NVOI
245 : SEGINI MLRCOE
246 : MLECOE.LECT(INOEU)=MLRCOE
247 : C
248 : C******* Loop involving the neighbors
249 : C We create the matrix to "pseudoinvert"
250 : C
251 : DO IVOI=1,NVOI,1
252 : NGV=MLENTI.LECT(IVOI)
253 : NLV=MLECEN.LECT(NGV)
254 : NLV1=MLELI1.LECT(NGV)
255 : IF((NLV .NE. 0) .OR. (NLV1 .NE. 0))THEN
256 : C
257 : C********** NGV is a 'centre' point or a Dirichlet point
258 : C
259 : IPCOOR=((NGV-1)*(IDIM+1))+1
260 : XV=MCOORD.XCOOR(IPCOOR)
261 : YV=MCOORD.XCOOR(IPCOOR+1)
262 : MATA.MAT(IVOI,1)=1
263 : MATA.MAT(IVOI,2)=XV-XS
264 : MATA.MAT(IVOI,3)=YV-YS
265 : IF(IDIM.EQ.3)THEN
266 : ZV=MCOORD.XCOOR(IPCOOR+2)
267 : MATA.MAT(IVOI,4)=ZV-ZS
268 : ENDIF
269 : MLRB.PROG(IVOI)=1.0D0
270 : ELSE
271 : NLV=MLELI2.LECT(NGV)
272 : IF(NLV .EQ. 0)THEN
273 : C
274 : C**************** NO B.C. specified on this point
275 : C
276 : WRITE(IOIMP,*) 'Boundary conditions ???'
276 : CALL ERREUR(21)
277 : GOTO 9999
278 : ELSE
279 : IPCOOR=((NGV-1)*(IDIM+1))+1
280 : XV=MCOORD.XCOOR(IPCOOR)
281 : YV=MCOORD.XCOOR(IPCOOR+1)
282 : DX=((XV - XS)**2)+((YV - YS)**2)
283 : IF(IDIM .EQ. 3)THEN
284 : ZV=MCOORD.XCOOR(IPCOOR+2)
285 : DX=DX+((ZV - ZS)**2)
286 : ENDIF
287 : NLN=MLEFAC.LECT(NGV)
288 : DX=DX**0.5D0
289 : MATA.MAT(IVOI,1)=0
290 : MATA.MAT(IVOI,2)=DX*MPNORM.VPOCHA(NLN,1)
291 : MATA.MAT(IVOI,3)=DX*MPNORM.VPOCHA(NLN,2)
292 : IF(IDIM .EQ. 3)THEN
293 : MATA.MAT(IVOI,4)=DX*MPNORM.VPOCHA(NLN,3)
294 : ENDIF
295 : MLRB.PROG(IVOI)=DX
296 : ENDIF
297 : ENDIF
298 : ENDDO
299 : C
300 : C TEST
301 : C do ivoi=1,nvoi,1
302 : C write(*,*)
303 : C & 'mata.mat(',ivoi,')=',(mata.mat(ivoi,i1),i1=1,nvmin,1)
304 : C write(*,*) 'b(',ivoi,')=',mlrb.prog(ivoi)
305 : C enddo
306 : C
307 : C
308 : C********** Now we have to invert this matrix
309 : C
310 : LOGSVD=.TRUE.
310 : CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT,
311 : & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD,
312 : & MLRTRA.PROG)
313 : IF(IERSVD.NE.0)THEN
314 : C
315 : C************* SVD decomposition of the matrix does not work
316 : C
317 : LOGSVD=.FALSE.
318 : ELSE
319 : C
320 : C************ We check the condition number of MATA
321 : C
322 : SMAX=0.0D0
323 : DO I1=1,NVMIN,1
324 : SMAX=MAX(SMAX,MLRSIG.PROG(I1))
325 : ENDDO
326 : SMIN=SMAX
327 : DO I1=1,NVMIN,1
328 : SMIN=MIN(SMIN,MLRSIG.PROG(I1))
329 : ENDDO
330 : IF((SMIN/SMAX) .LT. ERRTOL)THEN
331 : LOGSVD=.FALSE.
332 : ENDIF
333 : ENDIF
334 : C
335 : C TEST
336 : C write(*,*) 'LOGSVD=.FALSE.'
337 : C LOGSVD=.FALSE.
338 : C
339 : IF(LOGSVD)THEN
340 : C
341 : C********** INVSVD worked
342 : C
343 : C MATA = MATU MATSIG t(MATV)
344 : C inv(MATA) = MATV inv(MATSIG) t(MATU)
345 : C
346 : DO I4=1,NVMIN,1
347 : DO IVOI=1,NVOI,1
348 : DO I2=1,1,1
349 : C I2=1 is the only coefficient we are interested in
350 : MLRCOE.PROG(IVOI)=MLRCOE.PROG(IVOI)+
351 : & (MATV.MAT(I2,I4)*MATU.MAT(IVOI,I4)
352 : & *MLRB.PROG(IVOI)/MLRSIG.PROG(I4))
353 : ENDDO
354 : ENDDO
355 : ENDDO
356 : ELSE
357 : WRITE (IOIMP,*) 'rlenco.eso'
358 : C 22 0
359 : C Opération malvenue. Résultat douteux
359 : CALL ERREUR(22)
360 : C
361 : C********** INVSVD does not worked
362 : C For each neighbor k we have to compute the solution
363 : C of
364 : C
365 : C t(MATA) MATA x = t(MATA) * b
366 : C
367 : C where b= \sum_l e_l \delta(k,l) = e_k
368 : C
369 : C To do that, we compute
370 : C
371 : C X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
372 : C
373 : C X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
374 : C [t(MATA) MATA] X_0]
375 : C
376 : C
377 : C********** We compute [t(MATA) MATA]
378 : C We store it in the upper triangle of MTAA(NVMIN,NVMIN)
379 : C
380 : DO I1=1,NVMIN,1
381 : DO I2=I1,NVMIN,1
382 : MTAA.MAT(I1,I2)=0.0D0
383 : ENDDO
384 : ENDDO
385 :
386 : DO I1=1,NVMIN,1
387 : DO I2=I1,NVMIN,1
388 : DO I3=1,NVOI,1
389 : MTAA.MAT(I1,I2)=MTAA.MAT(I1,I2)+
390 : & (MATA.MAT(I3,I1)*MATA.MAT(I3,I2))
391 : ENDDO
392 : ENDDO
393 : ENDDO
394 : C
395 : C********** We compute [inve(t(MATA) MATA)]
396 : C CHOLDC stores it in the upper trianle of MINTAA(NVMIN,NVMIN)
397 : C
397 : CALL CHOLDC(NVMIN,NVMIN,MTAA.MAT,MLRIN1.PROG,MINTAA.MAT,
398 : & MLRIN2.PROG,IERR0)
399 : IF(IERR0.NE.0)THEN
400 : WRITE(IOIMP,*) 'subroutine rlenco.eso.'
401 : C 26 2
402 : C Tache impossible. Probablement données erronées
402 : CALL ERREUR(26)
403 : GOTO 9999
404 : ENDIF
405 : C
406 : C********** We complete MTAA and MINTAA
407 : C
408 : DO I1=1,NVMIN,1
409 : DO I2=I1+1,NVMIN,1
410 : MINTAA.MAT(I2,I1)=MINTAA.MAT(I1,I2)
411 : MTAA.MAT(I2,I1)=MTAA.MAT(I1,I2)
412 : ENDDO
413 : ENDDO
414 : C
415 : DO IVOI=1,NVOI,1
416 : C
417 : C************* We compute [t(MATA) . b] and we store it in MLRIN1.PROG
418 : C
419 : DO I1=1,NVMIN,1
420 : MLRIN1.PROG(I1)=MATA.MAT(IVOI,I1)*MLRB.PROG(IVOI)
421 : MLRIN2.PROG(I1)=0.0D0
422 : MLRIN3.PROG(I1)=0.0D0
423 : ENDDO
424 : C
425 : C************* X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
426 : C X_0(i1) into MLRIN2.PROG(I1)
427 : C
428 : DO I2=1,NVMIN,1
429 : DO I1=1,NVMIN,1
430 : MLRIN2.PROG(I1)=MLRIN2.PROG(I1)+
431 : & (MINTAA.MAT(I1,I2)*MLRIN1.PROG(I2))
432 : ENDDO
433 : ENDDO
434 : C
435 : C************ X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
436 : C [t(MATA) MATA] X_0]
437 : C
438 : C [t(MATA) MATA] X_0 into MLRIN3.PROG
439 : C
440 : DO I2=1,NVMIN,1
441 : DO I1=1,NVMIN,1
442 : MLRIN3.PROG(I1)=MLRIN3.PROG(I1)+
443 : & (MTAA.MAT(I1,I2)*MLRIN2.PROG(I2))
444 : ENDDO
445 : ENDDO
446 : C
447 : C************* Now we have
448 : C [t(MATA) . b] in MLRIN1.PROG
449 : C X_0(i1) in MLRIN2.PROG(I1)
450 : C [t(MATA) MATA] X_0 in MLRIN3.PROG
451 : C
452 : C X_1(i1) in MLRCOE.MAT(i1,IVOI)
453 : C
454 : DO I1=1,1,1
455 : C The only interesting one is I1=1
456 : DO I2=1,NVMIN,1
457 : MLRCOE.PROG(IVOI)=MLRCOE.PROG(IVOI)+
458 : & (MINTAA.MAT(I1,I2)*
459 : & (MLRIN1.PROG(I2)-MLRIN3.PROG(I2)))
460 : ENDDO
461 : MLRCOE.PROG(IVOI)=MLRCOE.PROG(IVOI)+
462 : & MLRIN2.PROG(I1)
463 : ENDDO
464 : ENDDO
465 : ENDIF
466 : CC
467 : CC TEST
468 : C write(*,*) 'ngs =', NGS
469 : C write(*,*) 'invide',LOGSVD
470 : C write(*,*) 'nvois =',(mlenti.lect(ivoi),ivoi=1,nvoi,1)
471 : C write(*,*) 'coeff =',(mlrcoe.prog(ivoi),ivoi=1,nvoi,1)
472 : C dx=0.0D0
473 : C do ivoi=1,nvoi,1
474 : C dx=dx+mlrcoe.prog(ivoi)
475 : C enddo
476 : C write(*,*) 'sum=',dx
477 : MLENTI=MLEPOI.LECT(INOEU)
478 : SEGDES MLENTI
479 : MLRCOE=MLECOE.LECT(INOEU)
480 : SEGDES MLRCOE
481 : ENDDO
482 : C
483 : SEGDES MELCEN
484 : SEGSUP MLECEN
485 : C
486 : IF(MELLI1 .NE. 0) SEGDES MELLI1
487 : IF(MELLI2 .NE. 0) SEGDES MELLI2
488 : SEGSUP MLELI1
489 : SEGSUP MLELI2
490 : C
491 : SEGDES MELSOM
492 : C
493 : SEGDES MLECOE
494 : SEGDES MLEPOI
495 : C
496 : SEGSUP MATA
497 : SEGSUP MATU
498 : SEGSUP MATV
499 : SEGSUP MLRSIG
500 : SEGSUP MLRTRA
501 : SEGSUP MLRB
502 : C
503 : SEGSUP MTAA
504 : SEGSUP MINTAA
505 : SEGSUP MLRIN1
506 : SEGSUP MLRIN2
507 : SEGSUP MLRIN3
508 : C
509 : SEGDES MPNORM
510 : SEGDES MELFAC
511 : SEGSUP MLEFAC
512 : C
513 : 9999 CONTINUE
514 : RETURN
515 : END
516 :
517 :
518 :
© Cast3M 2003 - All rights reserved.
Disclaimer