The TELEMAC-MASCARET system  trunk
om1101.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE om1101
3 ! *****************
4 !
5  &(op , dm,typdim,xm,typexm, dn,typdin,xn,typexn, c,
6  & nulone,nelbor,nbor,nelmax,ndiag,nptfr,nelebx,neleb)
7 !
8 !***********************************************************************
9 ! BIEF V7P0 21/08/2010
10 !***********************************************************************
11 !
12 !brief OPERATIONS ON MATRICES.
13 !code
14 !+ M: P1 TRIANGLE
15 !+ N: BOUNDARY MATRIX
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 TRANPOSE OF N TO M
26 !
27 !code
28 !+ CONVENTION FOR THE STORAGE OF EXTRA-DIAGONAL TERMS:
29 !+
30 !+ XM( ,1) ----> M(1,2)
31 !+ XM( ,2) ----> M(1,3)
32 !+ XM( ,3) ----> M(2,3)
33 !+ XM( ,4) ----> M(2,1)
34 !+ XM( ,5) ----> M(3,1)
35 !+ XM( ,6) ----> M(3,2)
36 !
37 !history J-M HERVOUET (LNHE)
38 !+ 23/06/2008
39 !+ V5P9
40 !+
41 !
42 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
43 !+ 13/07/2010
44 !+ V6P0
45 !+ Translation of French comments within the FORTRAN sources into
46 !+ English comments
47 !
48 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
49 !+ 21/08/2010
50 !+ V6P0
51 !+ Creation of DOXYGEN tags for automated documentation and
52 !+ cross-referencing of the FORTRAN sources
53 !
54 !history J-M HERVOUET (EDF LAB, LNHE)
55 !+ 13/03/2014
56 !+ V7P0
57 !+ Now written to enable different numbering of boundary points and
58 !+ boundary segments.
59 !
60 !history S.E.BOURBAN (HRW)
61 !+ 21/03/2017
62 !+ V7P3
63 !+ Replacement of the DATA declarations by the PARAMETER associates
64 !
65 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
66 !| C |-->| A GIVEN CONSTANT USED IN OPERATION OP
67 !| DM |<->| DIAGONAL OF M
68 !| DN |-->| DIAGONAL OF N
69 !| NBOR |-->| GLOBAL NUMBER OF BOUNDARY POINTS
70 !| NDIAG |-->| NUMBER OF TERMS IN THE DIAGONAL
71 !| NELBOR |-->| FOR THE KTH BOUNDARY EDGE, GIVES THE CORRESPONDING
72 !| | | ELEMENT.
73 !| NELEB |-->| NUMBER OF BOUNDARY ELEMENTS
74 !| NELEBX |-->| MAXIMUM NUMBER OF BOUNDARY ELEMENTS
75 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
76 !| NPTFR |-->| NUMBER OF BOUNDARY POINTS
77 !| NULONE |-->| GOES WITH ARRAY NELBOR. NELBOR GIVES THE
78 !| | | ADJACENT ELEMENT, NULONE GIVES THE LOCAL
79 !| | | NUMBER OF THE FIRST NODE OF THE BOUNDARY EDGE
80 !| | | I.E. 1, 2 OR 3 FOR TRIANGLES.
81 !| OP |-->| OPERATION TO BE DONE (SEE ABOVE)
82 !| TYPDIM |<->| TYPE OF DIAGONAL OF M:
83 !| | | TYPDIM = 'Q' : ANY VALUE
84 !| | | TYPDIM = 'I' : IDENTITY
85 !| | | TYPDIM = '0' : ZERO
86 !| TYPDIN |<->| TYPE OF DIAGONAL OF N:
87 !| | | TYPDIN = 'Q' : ANY VALUE
88 !| | | TYPDIN = 'I' : IDENTITY
89 !| | | TYPDIN = '0' : ZERO
90 !| TYPEXM |-->| TYPE OF OFF-DIAGONAL TERMS OF M:
91 !| | | TYPEXM = 'Q' : ANY VALUE
92 !| | | TYPEXM = 'S' : SYMMETRIC
93 !| | | TYPEXM = '0' : ZERO
94 !| TYPEXN |-->| TYPE OF OFF-DIAGONAL TERMS OF N:
95 !| | | TYPEXN = 'Q' : ANY VALUE
96 !| | | TYPEXN = 'S' : SYMMETRIC
97 !| | | TYPEXN = '0' : ZERO
98 !| XM |-->| OFF-DIAGONAL TERMS OF M
99 !| XN |-->| OFF-DIAGONAL TERMS OF N
100 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
101 !
102  USE bief, ex_om1101 => om1101
103 !
105  IMPLICIT NONE
106 !
107 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
108 !
109  INTEGER, INTENT(IN) :: NELMAX,NDIAG,NPTFR,NELEBX,NELEB
110  CHARACTER(LEN=8), INTENT(IN) :: OP
111  INTEGER, INTENT(IN) :: NBOR(*)
112  INTEGER, INTENT(IN) :: NULONE(nelebx),NELBOR(nelebx)
113  DOUBLE PRECISION, INTENT(IN) :: DN(*),XN(*)
114  DOUBLE PRECISION, INTENT(INOUT) :: DM(*),XM(nelmax,*)
115  CHARACTER(LEN=1), INTENT(INOUT) :: TYPDIM,TYPEXM,TYPDIN,TYPEXN
116  DOUBLE PRECISION, INTENT(IN) :: C
117 !
118 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
119 !
120  INTEGER K,IEL
121 !
122  DOUBLE PRECISION Z(1)
123 !
124 !-----------------------------------------------------------------------
125 !
126  INTEGER :: CORNSY(3,2)
127  parameter( cornsy = reshape( (/
128  & 1,3,5, 4,6,2 /), shape=(/ 3,2 /) ) )
129  INTEGER :: CORSYM(3)
130  parameter( corsym = (/ 1,3,2 /) )
131 !
132 !-----------------------------------------------------------------------
133 !
134  IF(op(1:8).EQ.'M=M+N ') THEN
135 !
136  IF(typdim.EQ.'Q'.AND.typdim.EQ.'Q'.AND.ndiag.GE.nptfr) THEN
137  CALL ovdb( 'X=X+Y ' , dm , dn , z , c , nbor , nptfr )
138  ELSE
139  WRITE(lu,199) typdim(1:1),op(1:8),typdin(1:1)
140 199 FORMAT(1x,'OM1101 (BIEF) : TYPDIM = ',a1,' NOT IMPLEMENTED',
141  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPDIN = ',a1)
142  CALL plante(0)
143  stop
144  ENDIF
145 !
146  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
147 !
148 ! CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
149 !
150  DO k = 1 , neleb
151  iel = nelbor(k)
152  xm( iel , cornsy(nulone(k),1) ) =
153  & xm( iel , cornsy(nulone(k),1) ) + xn(k)
154  xm( iel , cornsy(nulone(k),2) ) =
155  & xm( iel , cornsy(nulone(k),2) ) + xn(k+nelebx)
156  ENDDO
157 !
158  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
159 !
160 ! CASE WHERE M CAN BE ANYTHING AND N IS SYMMETRICAL
161 !
162  DO k = 1 , neleb
163  iel = nelbor(k)
164  xm( iel , cornsy(nulone(k),1) ) =
165  & xm( iel , cornsy(nulone(k),1) ) + xn(k)
166  xm( iel , cornsy(nulone(k),2) ) =
167  & xm( iel , cornsy(nulone(k),2) ) + xn(k)
168  ENDDO
169 !
170  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
171 !
172 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
173 !
174  DO k = 1 , neleb
175  iel = nelbor(k)
176  xm( iel , corsym(nulone(k)) ) =
177  & xm( iel , corsym(nulone(k)) ) + xn(k)
178  ENDDO
179 !
180  ELSE
181  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
182 99 FORMAT(1x,'OM1101 (BIEF) : TYPEXM = ',a1,' DOES NOT GO',
183  & /,1x,'FOR THE OPERATION : ',a8,' WITH TYPEXN = ',a1)
184  CALL plante(1)
185  stop
186  ENDIF
187 !
188 !-----------------------------------------------------------------------
189 !
190  ELSEIF(op(1:8).EQ.'M=M+TN ') THEN
191 !
192  CALL ovdb( 'X=X+Y ' , dm , dn , z , c , nbor , nptfr )
193 !
194  IF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'Q') THEN
195 !
196 ! CASE WHERE BOTH MATRICES ARE NONSYMMETRICAL
197 !
198  DO k = 1 , neleb
199  iel = nelbor(k)
200  xm( iel , cornsy(nulone(k),1) ) =
201  & xm( iel , cornsy(nulone(k),1) ) + xn(k+nelebx)
202  xm( iel , cornsy(nulone(k),2) ) =
203  & xm( iel , cornsy(nulone(k),2) ) + xn(k)
204  ENDDO
205 !
206  ELSEIF(typexm(1:1).EQ.'Q'.AND.typexn(1:1).EQ.'S') THEN
207 !
208 ! CASE WHERE N CAN BE ANYTHING AND N IS SYMMETRICAL
209 !
210  DO k = 1 , neleb
211  iel = nelbor(k)
212  xm( iel , cornsy(nulone(k),1) ) =
213  & xm( iel , cornsy(nulone(k),1) ) + xn(k)
214  xm( iel , cornsy(nulone(k),2) ) =
215  & xm( iel , cornsy(nulone(k),2) ) + xn(k)
216  ENDDO
217 !
218  ELSEIF(typexm(1:1).EQ.'S'.AND.typexn(1:1).EQ.'S') THEN
219 !
220 ! CASE WHERE BOTH MATRICES ARE SYMMETRICAL
221 !
222  DO k = 1 , neleb
223  iel = nelbor(k)
224  xm( iel , corsym(nulone(k)) ) =
225  & xm( iel , corsym(nulone(k)) ) + xn(k)
226  ENDDO
227 !
228  ELSE
229  WRITE(lu,99) typexm(1:1),op(1:8),typexn(1:1)
230  CALL plante(1)
231  stop
232  ENDIF
233 !
234 !-----------------------------------------------------------------------
235 !
236  ELSE
237 !
238  WRITE(lu,71) op
239 71 FORMAT(1x,'OM1101 (BIEF) : UNKNOWN OPERATION : ',a8)
240  CALL plante(1)
241  stop
242 !
243  ENDIF
244 !
245 !-----------------------------------------------------------------------
246 !
247  RETURN
248  END
subroutine ovdb(OP, X, Y, Z, C, NBOR, NPTFR)
Definition: ovdb.f:7
subroutine om1101(OP, DM, TYPDIM, XM, TYPEXM, DN, TYPDIN, XN, TYPEXN, C, NULONE, NELBOR, NBOR, NELMAX, NDIAG, NPTFR, NELEBX, NELEB)
Definition: om1101.f:8
Definition: bief.f:3