The TELEMAC-MASCARET system  trunk
omsegbor.f
Go to the documentation of this file.
1 ! *******************
2  SUBROUTINE omsegbor
3 ! *******************
4 !
5  &(op , dm,typdim,xm,typexm, dn,typdin,xn,typexn,c,
6  & ndiag,nseg1,nbor,nptfr,ielm1,ieln1,nseg11,
7  & iklbor,nelebx,neleb)
8 !
9 !***********************************************************************
10 ! BIEF V7P0 21/08/2010
11 !***********************************************************************
12 !
13 !brief OPERATIONS BETWEEN A MATRIX WITH EDGE-BASED STORAGE
14 !+ AND A BOUNDARY MATRIX.
15 !code
16 !+ D: DIAGONAL MATRIX
17 !+ C: CONSTANT
18 !+
19 !+ OP IS A STRING OF 8 CHARACTERS, WHICH INDICATES THE OPERATION TO BE
20 !+ PERFORMED ON MATRICES M AND N, D AND C.
21 !+
22 !+ THE RESULT IS MATRIX M.
23 !+
24 !+ OP = 'M=M+N ' : ADDS N TO M
25 !+ OP = 'M=M+TN ' : ADDS TRANSPOSE(N) TO M
26 !
27 !note IF BOTH MATRICES ARE QUADRATIC, THE NUMBER OF OFF-DIAGONAL TERMS
28 !+ IS MULTIPLIED BY 3 (THERE ARE 3 QUADRATIC SEGMENTS PER BOUNDARY
29 !+ SEGMENT), HENCE THE TERMS 3*NPTFR, WHICH ORIGINATES FROM THE FACT
30 !+ THAT SEGMENTS IN THE QUADRATIC TRIANGLE AND QUADRATIC SEGMENTS IN
31 !+ THE BOUNDARY SEGMENT ARE NUMBERED IN THE SAME ORDER.
32 !
33 !history J-M HERVOUET (LNHE)
34 !+ 12/02/2010
35 !+ V6P0
36 !+
37 !
38 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
39 !+ 13/07/2010
40 !+ V6P0
41 !+ Translation of French comments within the FORTRAN sources into
42 !+ English comments
43 !
44 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
45 !+ 21/08/2010
46 !+ V6P0
47 !+ Creation of DOXYGEN tags for automated documentation and
48 !+ cross-referencing of the FORTRAN sources
49 !
50 !history J-M HERVOUET (LNHE)
51 !+ 14/05/2012
52 !+ V6P2
53 !+ Bug corrected in the quadratic case.
54 !
55 !history J-M HERVOUET (EDF LAB, LNHE)
56 !+ 13/03/2014
57 !+ V7P0
58 !+ Now written to enable different numbering of boundary points and
59 !+ boundary segments.
60 !
61 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
62 !| C |-->| A GIVEN CONSTANT USED IN OPERATION OP
63 !| DM |<->| DIAGONAL OF M
64 !| DN |-->| DIAGONAL OF N
65 !| IELM1 |-->| TYPE OF ELEMENT OF M
66 !| IELN1 |-->| TYPE OF ELEMENT OF N
67 !| NBOR |-->| GLOBAL NUMBER OF BOUNDARY POINTS
68 !| NDIAG |-->| NUMBER OF TERMS IN THE DIAGONAL
69 !| NPTFR |-->| NUMBER OF BOUNDARY POINTS
70 !| NSEG1 |-->| NUMBER OF SEGMENTS CONSIDERED IN M
71 !| NSEG11 |-->| NUMBER OF LINEAR SEGMENTS
72 !| OP |-->| OPERATION TO BE DONE (SEE ABOVE)
73 !| TYPDIM |<->| TYPE OF DIAGONAL OF M:
74 !| | | TYPDIM = 'Q' : ANY VALUE
75 !| | | TYPDIM = 'I' : IDENTITY
76 !| | | TYPDIM = '0' : ZERO
77 !| TYPDIN |<->| TYPE OF DIAGONAL OF N:
78 !| | | TYPDIN = 'Q' : ANY VALUE
79 !| | | TYPDIN = 'I' : IDENTITY
80 !| | | TYPDIN = '0' : ZERO
81 !| TYPEXM |-->| TYPE OF OFF-DIAGONAL TERMS OF M:
82 !| | | TYPEXM = 'Q' : ANY VALUE
83 !| | | TYPEXM = 'S' : SYMMETRIC
84 !| | | TYPEXM = '0' : ZERO
85 !| TYPEXN |-->| TYPE OF OFF-DIAGONAL TERMS OF N:
86 !| | | TYPEXN = 'Q' : ANY VALUE
87 !| | | TYPEXN = 'S' : SYMMETRIC
88 !| | | TYPEXN = '0' : ZERO
89 !| XM |-->| OFF-DIAGONAL TERMS OF M
90 !| XN |-->| OFF-DIAGONAL TERMS OF N
91 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
92 !
93  USE bief, ex_omsegbor => omsegbor
94 !
95  USE declarations_special
96  IMPLICIT NONE
97 !
98 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
99 !
100  INTEGER, INTENT(IN) :: NDIAG,NSEG1,NPTFR,IELM1,IELN1,NSEG11
101  INTEGER, INTENT(IN) :: NELEBX,NELEB
102  INTEGER, INTENT(IN) :: NBOR(*),IKLBOR(nelebx,*)
103  CHARACTER(LEN=8), INTENT(IN) :: OP
104  DOUBLE PRECISION, INTENT(IN) :: DN(*),XN(nelebx,*)
105  DOUBLE PRECISION, INTENT(INOUT) :: DM(*),XM(nseg1,*)
106  CHARACTER(LEN=1), INTENT(INOUT) :: TYPDIM,TYPEXM,TYPDIN,TYPEXN
107  DOUBLE PRECISION, INTENT(IN) :: C
108 !
109 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
110 !
111  INTEGER IELEB,K
112  DOUBLE PRECISION Z(1)
113 !
114 !-----------------------------------------------------------------------
115 !
116  IF(op(1:8).EQ.'M=M+N ') THEN
117 !
118  IF(typdim.EQ.'Q'.AND.typdim.EQ.'Q'.AND.ndiag.GE.nptfr) THEN
119  CALL ovdb( 'X=X+Y ' , dm , dn , z , c , nbor , nptfr )
120 ! QUADRATIC POINTS IN THE MIDDLE OF SEGMENTS
121  IF(ielm1.EQ.13.AND.ieln1.EQ.2) THEN
122  DO ieleb=1,neleb
123  k=iklbor(ieleb,3)
124  dm(nbor(k))=dm(nbor(k))+dn(k)
125  ENDDO
126  ENDIF
127  ELSE
128  WRITE(lu,199) typdim(1:1),op(1:8),typdin(1:1)
129 199 FORMAT(1x,'OMSEGBOR (BIEF) : TYPDIM = ',a1,' NOT IMPLEMENTED',
130  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPDIN = ',a1)
131  CALL plante(1)
132  stop
133  ENDIF
134 !
135 ! THE BOUNDARY SEGMENTS ARE NUMBERED LIKE THE BOUNDARY NUMBERING
136 ! OF THEIR FIRST POINT (SEE STOSEG). HENCE THE (RELATIVELY SIMPLE)
137 ! IMPLEMENTATION BELOW. FURTHERMORE, ORISEG IS ALWAYS 1 FOR
138 ! BOUNDARY SEGMENTS, WHICH ALLOWS THE SHIFT OF NSEG11 AND 2*NSEG11
139 ! TO GET THE FIRST THEN THE SECOND HALF SEGMENT (SEE COMP_SEG).
140 !
141  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
142 !
143 ! CASE WHERE BOTH MATRICES ARE NON SYMMETRICAL
144  IF(ielm1.EQ.13.AND.ieln1.EQ.2) THEN
145 ! HERE XN(NELEBX,6) IN SEGMENTS POINT 3 IS THE MIDDLE
146 ! STORING IN XN : 1-2 1-3 2-3 2-1 3-1 2-3
147  CALL ov('X=X+Y ', x=xm( 1,1), y=xn(1,1),
148  & dim1=neleb)
149  CALL ov('X=X+Y ', x=xm( 1,2), y=xn(1,4),
150  & dim1=neleb)
151  CALL ov('X=X+Y ', x=xm( nseg11+1,1), y=xn(1,2),
152  & dim1=neleb)
153  CALL ov('X=X+Y ', x=xm( nseg11+1,2), y=xn(1,5),
154  & dim1=neleb)
155  CALL ov('X=X+Y ', x=xm(2*nseg11+1,1), y=xn(1,3),
156  & dim1=neleb)
157  CALL ov('X=X+Y ', x=xm(2*nseg11+1,2), y=xn(1,6),
158  & dim1=neleb)
159  ELSE
160 ! HERE XN(NELEBX,2)
161  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,1), dim1=neleb)
162  CALL ov('X=X+Y ', x=xm(1,2), y=xn(1,2), dim1=neleb)
163  ENDIF
164 !
165  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
166 !
167 ! CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
168  IF(ielm1.EQ.13.AND.ieln1.EQ.2) THEN
169 ! HERE XN(NELEBX,3)
170  CALL ov('X=X+Y ', x=xm( 1,1), y=xn(1,1),
171  & dim1=neleb)
172  CALL ov('X=X+Y ', x=xm( 1,2), y=xn(1,1),
173  & dim1=neleb)
174  CALL ov('X=X+Y ', x=xm( nseg11+1,1), y=xn(1,2),
175  & dim1=neleb)
176  CALL ov('X=X+Y ', x=xm( nseg11+1,2), y=xn(1,2),
177  & dim1=neleb)
178  CALL ov('X=X+Y ', x=xm(2*nseg11+1,1), y=xn(1,3),
179  & dim1=neleb)
180  CALL ov('X=X+Y ', x=xm(2*nseg11+1,2), y=xn(1,3),
181  & dim1=neleb)
182  ELSE
183 ! HERE XN(NPTFR,1)
184  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,1), dim1=neleb)
185  CALL ov('X=X+Y ', x=xm(1,2), y=xn(1,1), dim1=neleb)
186  ENDIF
187 !
188  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
189 !
190 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
191  IF(ielm1.EQ.13.AND.ieln1.EQ.2) THEN
192 ! HERE XN(NELEBX,3)
193  CALL ov('X=X+Y ', x=xm( 1,1), y=xn(1,1),
194  & dim1=neleb)
195  CALL ov('X=X+Y ', x=xm( nseg11+1,1), y=xn(1,2),
196  & dim1=neleb)
197  CALL ov('X=X+Y ', x=xm(2*nseg11+1,1), y=xn(1,3),
198  & dim1=neleb)
199  ELSE
200 ! HERE XN(NELEBX,1)
201  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,1), dim1=neleb)
202  ENDIF
203 !
204  ELSE
205 !
206  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
207 99 FORMAT(1x,'OMSEGBOR (BIEF) : TYPEXM = ',a1,' DOES NOT GO',
208  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPEXN = ',a1)
209  CALL plante(1)
210  stop
211 !
212  ENDIF
213 !
214 !-----------------------------------------------------------------------
215 !
216  ELSEIF(op(1:8).EQ.'M=M+TN ') THEN
217 !
218  IF(typdim.EQ.'Q'.AND.typdim.EQ.'Q'.AND.ndiag.GE.nptfr) THEN
219  CALL ovdb( 'X=X+Y ' , dm , dn , z , c , nbor , nptfr )
220 ! QUADRATIC POINTS IN THE MIDDLE OF SEGMENTS
221  IF(ielm1.EQ.13.AND.ieln1.EQ.2) THEN
222  DO ieleb=1,neleb
223  k=iklbor(ieleb,3)
224  dm(nbor(k))=dm(nbor(k))+dn(k)
225  ENDDO
226  ENDIF
227  ELSE
228  WRITE(lu,199) typdim(1:1),op(1:8),typdin(1:1)
229  CALL plante(1)
230  stop
231  ENDIF
232 !
233  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
234 !
235 ! CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
236  IF(ielm1.EQ.13.AND.ieln1.EQ.2) THEN
237 ! HERE XN(NELEBX,6)
238  CALL ov('X=X+Y ', x=xm( 1,1), y=xn(1,4),
239  & dim1=neleb)
240  CALL ov('X=X+Y ', x=xm( 1,2), y=xn(1,1),
241  & dim1=neleb)
242  CALL ov('X=X+Y ', x=xm( nseg11+1,1), y=xn(1,5),
243  & dim1=neleb)
244  CALL ov('X=X+Y ', x=xm( nseg11+1,2), y=xn(1,2),
245  & dim1=neleb)
246  CALL ov('X=X+Y ', x=xm(2*nseg11+1,1), y=xn(1,6),
247  & dim1=neleb)
248  CALL ov('X=X+Y ', x=xm(2*nseg11+1,2), y=xn(1,3),
249  & dim1=neleb)
250  ELSE
251  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,2), dim1=neleb)
252  CALL ov('X=X+Y ', x=xm(1,2), y=xn(1,1), dim1=neleb)
253  ENDIF
254 !
255  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
256 !
257 ! CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
258  IF(ielm1.EQ.13.AND.ieln1.EQ.2) THEN
259 ! HERE XN(NELEBX,3)
260  CALL ov('X=X+Y ', x=xm( 1,1), y=xn(1,1),
261  & dim1=neleb)
262  CALL ov('X=X+Y ', x=xm( 1,2), y=xn(1,1),
263  & dim1=neleb)
264  CALL ov('X=X+Y ', x=xm( nseg11+1,1), y=xn(1,2),
265  & dim1=neleb)
266  CALL ov('X=X+Y ', x=xm( nseg11+1,2), y=xn(1,2),
267  & dim1=neleb)
268  CALL ov('X=X+Y ', x=xm(2*nseg11+1,1), y=xn(1,3),
269  & dim1=neleb)
270  CALL ov('X=X+Y ', x=xm(2*nseg11+1,2), y=xn(1,3),
271  & dim1=neleb)
272  ELSE
273  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,1), dim1=neleb)
274  CALL ov('X=X+Y ', x=xm(1,2), y=xn(1,1), dim1=neleb)
275  ENDIF
276 !
277  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
278 !
279 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
280  IF(ielm1.EQ.13.AND.ieln1.EQ.2) THEN
281 ! HERE XN(NELEBX,3)
282  CALL ov('X=X+Y ', x=xm( 1,1), y=xn(1,1),
283  & dim1=neleb)
284  CALL ov('X=X+Y ', x=xm( nseg11+1,1), y=xn(1,2),
285  & dim1=neleb)
286  CALL ov('X=X+Y ', x=xm(2*nseg11+1,1), y=xn(1,3),
287  & dim1=neleb)
288  ELSE
289  CALL ov('X=X+Y ', x=xm(1,1), y=xn(1,1), dim1=neleb)
290  ENDIF
291 !
292  ELSE
293 !
294  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
295  CALL plante(1)
296  stop
297 !
298  ENDIF
299 !
300 !-----------------------------------------------------------------------
301 !
302  ELSE
303 !
304  WRITE(lu,71) op
305 71 FORMAT(1x,'OMSEGBOR (BIEF) : UNKNOWN OPERATION : ',a8)
306  CALL plante(1)
307  stop
308 !
309  ENDIF
310 !
311 !-----------------------------------------------------------------------
312 !
313  RETURN
314  END
ov
subroutine ov(OP, X, Y, Z, C, DIM1)
Definition: ov.f:7
declarations_special
Definition: declarations_special.F:3
omsegbor
subroutine omsegbor(OP, DM, TYPDIM, XM, TYPEXM, DN, TYPDIN, XN, TYPEXN, C, NDIAG, NSEG1, NBOR, NPTFR, IELM1, IELN1, NSEG11, IKLBOR, NELEBX, NELEB)
Definition: omsegbor.f:9
declarations_special::lu
integer lu
Definition: declarations_special.F:45
ovdb
subroutine ovdb(OP, X, Y, Z, C, NBOR, NPTFR)
Definition: ovdb.f:7
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:21 for The TELEMAC-MASCARET system by doxygen 1.8.13