The TELEMAC-MASCARET system  trunk
diri09.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE diri09
3 ! *****************
4 !
5  &(x1,x2,x3,
6  & a11,a12,a13,a21,a22,a23,a31,a32,a33,
7  & sm1,sm2,sm3,t1,t2,t3,t4,t5,t6,
8  & xbor1,xbor2,xbor3,lidir1,lidir2,lidir3,
9  & mesh,kdir,msk,maskpt)
10 !
11 !***********************************************************************
12 ! BIEF V6P1 21/08/2010
13 !***********************************************************************
14 !
15 !brief TREATS THE DIRICHLET POINTS FOR THE FOLLOWING
16 !+ SYSTEM (BLOCK OF 9 MATRICES):
17 !code
18 !+ ( A11 A12 A13 ) ( X1 ) ( SM1 )
19 !+ ( ) ( ) ( )
20 !+ ( T ) ( ) ( )
21 !+ ( A21 A22 A23 ) ( X2 ) = ( SM2 )
22 !+ ( ) ( ) ( )
23 !+ ( T T ) ( ) ( )
24 !+ ( A31 A32 A33 ) ( X3 ) ( SM3 )
25 !
26 !note TRANSPOSING A21 A31 AND A32 MAKES IT POSSIBLE TO USE ONLY
27 !+ ONE CALL FOR A12 AND A21, A31 AND A13, A32 AND A23 WHEN
28 !+ THE BLOCK IS SYMMETRICAL.
29 !
30 !history J-M HERVOUET (LNH)
31 !+ 30/01/95
32 !+ V5P1
33 !+
34 !
35 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
36 !+ 13/07/2010
37 !+ V6P0
38 !+ Translation of French comments within the FORTRAN sources into
39 !+ English comments
40 !
41 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
42 !+ 21/08/2010
43 !+ V6P0
44 !+ Creation of DOXYGEN tags for automated documentation and
45 !+ cross-referencing of the FORTRAN sources
46 !
47 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
48 !| A12 |<->| MATRIX IN THE 3x3 LINEAR SYSTEM
49 !| A13 |<->| MATRIX IN THE 3x3 LINEAR SYSTEM
50 !| A21 |<->| MATRIX IN THE 3x3 LINEAR SYSTEM
51 !| A22 |<->| MATRIX IN THE 3x3 LINEAR SYSTEM
52 !| A23 |<->| MATRIX IN THE 3x3 LINEAR SYSTEM
53 !| A31 |<->| MATRIX IN THE 3x3 LINEAR SYSTEM
54 !| A32 |<->| MATRIX IN THE 3x3 LINEAR SYSTEM
55 !| A33 |<->| MATRIX IN THE 3x3 LINEAR SYSTEM
56 !| KDIR |-->| CONVENTION FOR DIRICHLET BOUNDARY CONDITIONS
57 !| LIDIR1 |-->| TYPES OF BOUNDARY CONDITIONS FOR VARIABLE 1
58 !| | | IF LIMDIR(K) = KDIR LE KTH BOUNDARY POINT
59 !| | | IS OF DIRICHLET TYPE.
60 !| LIDIR2 |-->| TYPES OF BOUNDARY CONDITIONS FOR VARIABLE 2
61 !| | | IF LIMDIR(K) = KDIR THE KTH BOUNDARY POINT
62 !| | | IS OF DIRICHLET TYPE.
63 !| LIDIR2 |-->| TYPES OF BOUNDARY CONDITIONS FOR VARIABLE 2
64 !| | | IF LIMDIR(K) = KDIR THE KTH BOUNDARY POINT
65 !| | | IS OF DIRICHLET TYPE.
66 !| LIDIR3 |-->| TYPES OF BOUNDARY CONDITIONS FOR VARIABLE 3
67 !| | | IF LIMDIR(K) = KDIR THE KTH BOUNDARY POINT
68 !| | | IS OF DIRICHLET TYPE.
69 !| MASKPT |-->| MASKING PER POINT.
70 !| | | =1. : NORMAL =0. : MASKED
71 !| MESH |-->| MESH STRUCTURE
72 !| MSK |-->| IF YES, THERE IS MASKED ELEMENTS.
73 !| SM1 |-->| FIRST RIGHT-HAND SIDE OF THE SYSTEM.
74 !| SM2 |-->| SECOND RIGHT-HAND SIDE OF THE SYSTEM.
75 !| SM3 |-->| THIRD RIGHT-HAND SIDE OF THE SYSTEM.
76 !| T1 |<->| WORK DOUBLE PRECISION ARRAY IN A BIEF_OBJ
77 !| T2 |<->| WORK DOUBLE PRECISION ARRAY IN A BIEF_OBJ
78 !| T3 |<->| WORK DOUBLE PRECISION ARRAY IN A BIEF_OBJ
79 !| T4 |<->| WORK DOUBLE PRECISION ARRAY IN A BIEF_OBJ
80 !| T5 |<->| WORK DOUBLE PRECISION ARRAY IN A BIEF_OBJ
81 !| T6 |<->| WORK DOUBLE PRECISION ARRAY IN A BIEF_OBJ
82 !| XBOR1 |-->| DIRICHLET BOUNDARY CONDITIONS OF VARIABLE 1
83 !| XBOR2 |-->| DIRICHLET BOUNDARY CONDITIONS OF VARIABLE 2
84 !| XBOR3 |-->| DIRICHLET BOUNDARY CONDITIONS OF VARIABLE 3
85 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
86 !
87  USE bief, ex_diri09 => diri09
88 !
90  IMPLICIT NONE
91 !
92 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
93 !
94  TYPE(bief_obj), INTENT(INOUT) :: X1,X2,X3,SM1,SM2,SM3
95  TYPE(bief_obj), INTENT(INOUT) :: T1,T2,T3,T4,T5,T6
96  TYPE(bief_obj), INTENT(INOUT) :: A11,A12,A13,A21,A22
97  TYPE(bief_obj), INTENT(INOUT) :: A23,A31,A32,A33
98  TYPE(bief_obj), INTENT(IN) :: XBOR1,XBOR2,XBOR3,MASKPT
99  INTEGER, INTENT(IN) :: LIDIR1(*),LIDIR2(*),LIDIR3(*)
100  INTEGER, INTENT(IN) :: KDIR
101  TYPE(bief_mesh), INTENT(INOUT):: MESH
102  LOGICAL, INTENT(IN) :: MSK
103 !
104 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
105 !
106  DOUBLE PRECISION C
107 !
108  CHARACTER(LEN=1) STODIA
109 !
110 !-----------------------------------------------------------------------
111 !
112 ! 1) BUILDS ARRAYS T1,T2,T3 CONTAINING:
113 ! THE X1, X2 AND X3 IMPOSED VALUES IF THE POINT IS OF TYPE DIRICHLET
114 ! 0 OTHERWISE
115 !
116 ! X1,X2,X3 TAKE THEIR DIRICHLET VALUE
117 !
118 !=======================================================================
119 !
120 ! BOUNDARY CONDITION FOR X1 : "XBOR1" IMPOSED
121 !
122  CALL cpstvc(x1,t1)
123  CALL os('X=C ', x=t1, c=0.d0)
124  CALL osdbif ( 'X=Y ',t1,xbor1,lidir1,kdir,mesh)
125 !
126 !-----------------------------------------------------------------------
127 !
128 ! BOUNDARY CONDITION FOR X2 : "XBOR2" IMPOSED
129 !
130  CALL cpstvc(x2,t2)
131  CALL os('X=C ', x=t2, c=0.d0)
132  CALL osdbif ( 'X=Y ',t2,xbor2,lidir2,kdir,mesh)
133 !
134 !-----------------------------------------------------------------------
135 !
136 ! BOUNDARY CONDITION FOR X3 : "XBOR3" IMPOSED
137 !
138  CALL cpstvc(x3,t3)
139  CALL os('X=C ', x=t3, c=0.d0)
140  CALL osdbif ( 'X=Y ',t3,xbor3,lidir3,kdir,mesh)
141 !
142 !=======================================================================
143 !
144 ! 2) COMPUTES THE PRODUCT OF THE MATRIX FOR THE SYSTEM TO SOLVE
145 ! AND T1,T2,T3
146 ! THE RESULT IS DEDUCTED FROM THE SECOND MEMBERS
147 !
148  CALL matvec('X=AY ',t4,a11,t1,c,mesh,lego=.false.)
149  CALL matvec('X=X+AY ',t4,a12,t2,c,mesh,lego=.false.)
150  CALL matvec('X=X+AY ',t4,a13,t3,c,mesh,lego=.true. )
151  CALL matvec('X=AY ',t5,a21,t1,c,mesh,lego=.false.)
152  CALL matvec('X=X+AY ',t5,a22,t2,c,mesh,lego=.false.)
153  CALL matvec('X=X+AY ',t5,a23,t3,c,mesh,lego=.true. )
154  CALL matvec('X=AY ',t6,a31,t1,c,mesh,lego=.false.)
155  CALL matvec('X=X+AY ',t6,a32,t2,c,mesh,lego=.false.)
156  CALL matvec('X=X+AY ',t6,a33,t3,c,mesh,lego=.true. )
157 !
158  CALL cpstvc(x1,sm1)
159  CALL cpstvc(x2,sm2)
160  CALL cpstvc(x3,sm3)
161  CALL os('X=X-Y ', x=sm1, y=t4)
162  CALL os('X=X-Y ', x=sm2, y=t5)
163  CALL os('X=X-Y ', x=sm3, y=t6)
164 !
165 !=======================================================================
166 !
167 ! SECOND MEMBERS OF THE EQUATIONS FOR DIRICHLET POINTS
168 ! PREPARES THE LINEAR SYSTEM
169 !
170  CALL diraux(sm1,a11%D,xbor1,t1,x1,lidir1,kdir,mesh)
171  CALL diraux(sm2,a22%D,xbor2,t2,x2,lidir2,kdir,mesh)
172  CALL diraux(sm3,a33%D,xbor3,t3,x3,lidir3,kdir,mesh)
173 !
174 ! CALLS OV RATHER THAN OS BECAUSE SM1 AND MASKPT DON'T ALWAYS
175 ! HAVE THE SAME LENGTH
176 !
177  IF(msk) THEN
178  CALL ov('X=XY ', x=sm1%R, y=maskpt%R, dim1=sm1%DIM1)
179  CALL ov('X=XY ', x=x1%R, y=maskpt%R, dim1=x1%DIM1)
180  CALL ov('X=XY ', x=t1%R, y=maskpt%R, dim1=t1%DIM1)
181  CALL ov('X=XY ', x=sm2%R, y=maskpt%R, dim1=sm2%DIM1)
182  CALL ov('X=XY ', x=x2%R, y=maskpt%R, dim1=x2%DIM1)
183  CALL ov('X=XY ', x=t2%R, y=maskpt%R, dim1=t2%DIM1)
184  CALL ov('X=XY ', x=sm3%R, y=maskpt%R, dim1=sm3%DIM1)
185  CALL ov('X=XY ', x=x3%R, y=maskpt%R, dim1=x3%DIM1)
186  CALL ov('X=XY ', x=t3%R, y=maskpt%R, dim1=t3%DIM1)
187  ENDIF
188 !
189 !=======================================================================
190 !
191 ! ERASES THE LINES AND COLUMNS FOR DIRICHLET POINTS
192 !
193 ! IT'S EQUIVALENT TO A DIAGONAL PRECONDITIONING WITH ARRAYS
194 ! T1,T2,T3
195 !
196 ! DOES NOT ALTER A11,A22,A33 DIAGONALS
197 ! BY GIVING THEM A DUMMY TYPE : '0'
198 !
199 !
200 !=======================================================================
201 ! A11 PRECONDITIONING :
202 !=======================================================================
203 !
204  stodia = a11%TYPDIA
205  a11%TYPDIA='0'
206  CALL om('M=DMD ', m=a11, d=t1, mesh=mesh)
207  a11%TYPDIA=stodia
208 !
209 !=======================================================================
210 ! A12 PRECONDITIONING :
211 !=======================================================================
212 !
213  CALL om('M=DM ', m=a12, d=t1, mesh=mesh)
214  CALL om('M=MD ', m=a12, d=t2, mesh=mesh)
215 !
216 !=======================================================================
217 ! A13 PRECONDITIONING :
218 !=======================================================================
219 !
220  CALL om('M=DM ', m=a13, d=t1, mesh=mesh)
221  CALL om('M=MD ', m=a13, d=t3, mesh=mesh)
222 !
223 !=======================================================================
224 ! A21 PRECONDITIONING :
225 !=======================================================================
226 !
227  CALL om('M=DM ', m=a21, d=t2, mesh=mesh)
228  CALL om('M=MD ', m=a21, d=t1, mesh=mesh)
229 !
230 !=======================================================================
231 ! A22 PRECONDITIONING :
232 !=======================================================================
233 !
234  stodia = a22%TYPDIA
235  a22%TYPDIA='0'
236  CALL om('M=DMD ', m=a22, d=t2, mesh=mesh)
237  a22%TYPDIA=stodia
238 !
239 !=======================================================================
240 ! A23 PRECONDITIONING :
241 !=======================================================================
242 !
243  CALL om('M=DM ', m=a23, d=t2, mesh=mesh)
244  CALL om('M=MD ', m=a23, d=t3, mesh=mesh)
245 !
246 !=======================================================================
247 ! A31 PRECONDITIONING :
248 !=======================================================================
249 !
250  CALL om('M=DM ', m=a31, d=t3, mesh=mesh)
251  CALL om('M=MD ', m=a31, d=t1, mesh=mesh)
252 !
253 !=======================================================================
254 ! A32 PRECONDITIONING :
255 !=======================================================================
256 !
257  CALL om('M=DM ', m=a32, d=t3, mesh=mesh)
258  CALL om('M=MD ', m=a32, d=t2, mesh=mesh)
259 !
260 !=======================================================================
261 ! A33 PRECONDITIONING :
262 !=======================================================================
263 !
264  stodia = a33%TYPDIA
265  a33%TYPDIA='0'
266  CALL om('M=DMD ', m=a33, d=t3, mesh=mesh)
267  a33%TYPDIA=stodia
268 !
269 !-----------------------------------------------------------------------
270 !
271  RETURN
272  END
subroutine ov(OP, X, Y, Z, C, DIM1)
Definition: ov.f:7
subroutine diraux(X, Y, Z, W, F, INDIC, CRITER, MESH)
Definition: diraux.f:7
subroutine diri09(X1, X2, X3, A11, A12, A13, A21, A22, A23, A31, A32, A33, SM1, SM2, SM3, T1, T2, T3, T4, T5, T6, XBOR1, XBOR2, XBOR3, LIDIR1, LIDIR2, LIDIR3, MESH, KDIR, MSK, MASKPT)
Definition: diri09.f:11
subroutine om(OP, M, N, D, C, MESH)
Definition: om.f:7
subroutine osdbif(OP, X, Y, INDIC, CRITER, MESH)
Definition: osdbif.f:7
subroutine cpstvc(X, Y)
Definition: cpstvc.f:7
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
subroutine matvec(OP, X, A, Y, C, MESH, LEGO)
Definition: matvec.f:7
Definition: bief.f:3