SFEMaNS  version 5.3
Reference documentation for SFEMaNS
assembling_rhs.f90
Go to the documentation of this file.
1 MODULE rhs_para_assembling
2 
3 CONTAINS
4 
5  SUBROUTINE rhs_3x3(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_tensor_scaln_bdy)
6  USE def_type_mesh
7  USE my_util
8 #include "petsc/finclude/petsc.h"
9  USE petsc
10  IMPLICIT NONE
11  TYPE(mesh_type) :: mesh
12  TYPE(petsc_csr_la) :: vv_3_LA
13  INTEGER, INTENT(IN) :: mode
14  REAL(KIND=8), DIMENSION(mesh%dom_me*mesh%gauss%l_G,6), INTENT(IN) :: rhs_in
15  REAL(KIND=8), DIMENSION(3,mesh%dom_me*mesh%gauss%l_G,6), INTENT(IN), OPTIONAL :: opt_tensor
16  REAL(KIND=8), DIMENSION(:,:), INTENT(IN), OPTIONAL :: opt_tensor_scaln_bdy
17  REAL(KIND=8), DIMENSION(3*mesh%gauss%n_w) :: v145_loc, v236_loc
18  INTEGER, DIMENSION(3*mesh%gauss%n_w) :: idxm
19  INTEGER, DIMENSION(3*mesh%gauss%n_ws) :: idxms
20  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
21  INTEGER, DIMENSION(mesh%gauss%n_ws) :: js_loc
22  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
23  REAL(KIND=8), DIMENSION(3,6) :: tensor_loc
24  REAL(KIND=8), DIMENSION(6) :: smb, f_tensor
25  INTEGER :: m, l, i, ni, index, nw, ix, ki, iglob, type, k, ms, ls, indexs, nws
26  REAL(KIND=8) :: ray
27  vec :: vb_145, vb_236
28  petscerrorcode :: ierr
29 
30  CALL vecset(vb_145, 0.d0, ierr)
31  CALL vecset(vb_236, 0.d0, ierr)
32 
33  nw = mesh%gauss%n_w
34  index = 0
35  DO m = 1, mesh%dom_me
36  j_loc = mesh%jj(:,m)
37 
38  DO ki = 1, 3
39  DO ni = 1, nw
40  i = j_loc(ni)
41  iglob = vv_3_la%loc_to_glob(ki,i)
42  ix = (ki-1)*nw+ni
43  idxm(ix) = iglob-1
44  END DO
45  END DO
46 
47  v145_loc = 0.d0
48  v236_loc = 0.d0
49  DO l = 1, mesh%gauss%l_G
50  index = index + 1
51  dw_loc = mesh%gauss%dw(:,:,l,m)
52 
53  !===Compute radius of Gauss point
54  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
55 
56  IF (PRESENT(opt_tensor)) THEN
57  DO TYPE = 1, 6
58  DO k = 1, 3
59  tensor_loc(k,type) = opt_tensor(k,index,type)
60  END DO
61  END DO
62 
63  f_tensor(1) = (-mode*tensor_loc(1,4) + tensor_loc(2,3))/ray
64  f_tensor(2) = (mode*tensor_loc(1,3) + tensor_loc(2,4))/ray
65  f_tensor(3) = (-tensor_loc(1,3) - mode*tensor_loc(2,4))/ray
66  f_tensor(4) = (-tensor_loc(1,4) + mode*tensor_loc(2,3))/ray
67  f_tensor(5) = -tensor_loc(3,4)*mode/ray
68  f_tensor(6) = tensor_loc(3,3)*mode/ray
69  ELSE
70  f_tensor = 0.d0
71  END IF
72 
73  smb = (rhs_in(index,:)+f_tensor)*ray*mesh%gauss%rj(l,m)
74 
75  ki = 1
76  DO ni = 1, nw
77  ix = (ki-1)*nw+ni
78  v145_loc(ix) = v145_loc(ix) + mesh%gauss%ww(ni,l)*smb(1)
79  v236_loc(ix) = v236_loc(ix) + mesh%gauss%ww(ni,l)*smb(2)
80  END DO
81 
82  ki = 2
83  DO ni = 1, nw
84  ix = (ki-1)*nw+ni
85  v145_loc(ix) = v145_loc(ix) + mesh%gauss%ww(ni,l)*smb(4)
86  v236_loc(ix) = v236_loc(ix) + mesh%gauss%ww(ni,l)*smb(3)
87  END DO
88 
89  ki = 3
90  DO ni = 1, nw
91  ix = (ki-1)*nw+ni
92  v145_loc(ix) = v145_loc(ix) + mesh%gauss%ww(ni,l)*smb(5)
93  v236_loc(ix) = v236_loc(ix) + mesh%gauss%ww(ni,l)*smb(6)
94  END DO
95 
96  IF (PRESENT(opt_tensor)) THEN
97  ki = 1
98  DO ni = 1, nw
99  ix = (ki-1)*nw+ni
100  v145_loc(ix) = v145_loc(ix) &
101  + (dw_loc(1,ni)*tensor_loc(1,1) + dw_loc(2,ni)*tensor_loc(1,5))*ray*mesh%gauss%rj(l,m)
102  v236_loc(ix) = v236_loc(ix) &
103  + (dw_loc(1,ni)*tensor_loc(1,2) + dw_loc(2,ni)*tensor_loc(1,6))*ray*mesh%gauss%rj(l,m)
104  END DO
105  ki = 2
106  DO ni = 1, nw
107  ix = (ki-1)*nw+ni
108  v145_loc(ix) = v145_loc(ix) &
109  + (dw_loc(1,ni)*tensor_loc(2,2) + dw_loc(2,ni)*tensor_loc(2,6))*ray*mesh%gauss%rj(l,m)
110  v236_loc(ix) = v236_loc(ix) &
111  + (dw_loc(1,ni)*tensor_loc(2,1) + dw_loc(2,ni)*tensor_loc(2,5))*ray*mesh%gauss%rj(l,m)
112  END DO
113  ki = 3
114  DO ni = 1, nw
115  ix = (ki-1)*nw+ni
116  v145_loc(ix) = v145_loc(ix) &
117  + (dw_loc(1,ni)*tensor_loc(3,1) + dw_loc(2,ni)*tensor_loc(3,5))*ray*mesh%gauss%rj(l,m)
118  v236_loc(ix) = v236_loc(ix) &
119  + (dw_loc(1,ni)*tensor_loc(3,2) + dw_loc(2,ni)*tensor_loc(3,6))*ray*mesh%gauss%rj(l,m)
120  END DO
121  END IF
122  ENDDO !l_G
123 
124  CALL vecsetvalues(vb_145, 3*nw, idxm, v145_loc, add_values, ierr)
125  CALL vecsetvalues(vb_236, 3*nw, idxm, v236_loc, add_values, ierr)
126 
127  ENDDO !dom_me
128 
129  IF (PRESENT(opt_tensor_scaln_bdy)) THEN
130  nws = mesh%gauss%n_ws
131  indexs = 0
132  DO ms = 1, mesh%dom_mes
133  js_loc = mesh%jjs(:,ms)
134 
135  DO ki = 1, 3
136  DO ni = 1, nws
137  i = js_loc(ni)
138  iglob = vv_3_la%loc_to_glob(ki,i)
139  ix = (ki-1)*nws+ni
140  idxms(ix) = iglob-1
141  END DO
142  END DO
143 
144  v145_loc = 0.d0
145  v236_loc = 0.d0
146  DO ls = 1, mesh%gauss%l_Gs
147  indexs = indexs + 1
148 
149  !===Compute radius of Gauss point
150  ray = sum(mesh%rr(1,js_loc)*mesh%gauss%wws(:,ls))
151 
152  !===Don't integrate on r=0
153  IF (ray.LT.1.d-10) cycle
154 
155  smb(:) = opt_tensor_scaln_bdy(indexs,:)*ray*mesh%gauss%rjs(ls,ms)
156 
157  ki = 1
158  DO ni = 1, nws
159  ix = (ki-1)*nws+ni
160  v145_loc(ix) = v145_loc(ix) + mesh%gauss%wws(ni,ls)*smb(1)
161  v236_loc(ix) = v236_loc(ix) + mesh%gauss%wws(ni,ls)*smb(2)
162  END DO
163 
164  ki = 2
165  DO ni = 1, nws
166  ix = (ki-1)*nws+ni
167  v145_loc(ix) = v145_loc(ix) + mesh%gauss%wws(ni,ls)*smb(4)
168  v236_loc(ix) = v236_loc(ix) + mesh%gauss%wws(ni,ls)*smb(3)
169  END DO
170 
171  ki = 3
172  DO ni = 1, nws
173  ix = (ki-1)*nws+ni
174  v145_loc(ix) = v145_loc(ix) + mesh%gauss%wws(ni,ls)*smb(5)
175  v236_loc(ix) = v236_loc(ix) + mesh%gauss%wws(ni,ls)*smb(6)
176  END DO
177 
178 
179  END DO
180  CALL vecsetvalues(vb_145, 3*nws, idxms, v145_loc, add_values, ierr)
181  CALL vecsetvalues(vb_236, 3*nws, idxms, v236_loc, add_values, ierr)
182  END DO
183  END IF
184 
185  CALL vecassemblybegin(vb_145,ierr)
186  CALL vecassemblyend(vb_145,ierr)
187  CALL vecassemblybegin(vb_236,ierr)
188  CALL vecassemblyend(vb_236,ierr)
189  END SUBROUTINE rhs_3x3
190 
191 
192  SUBROUTINE rhs_3x3_art_comp(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_grad_div)
193  USE def_type_mesh
194  USE my_util
195 #include "petsc/finclude/petsc.h"
196  USE petsc
197  IMPLICIT NONE
198  TYPE(mesh_type) :: mesh
199  TYPE(petsc_csr_la) :: vv_3_LA
200  INTEGER, INTENT(IN) :: mode
201  REAL(KIND=8), DIMENSION(mesh%dom_me*mesh%gauss%l_G,6), INTENT(IN) :: rhs_in
202  REAL(KIND=8), DIMENSION(3,mesh%dom_me*mesh%gauss%l_G,6), INTENT(IN), OPTIONAL :: opt_tensor
203  REAL(KIND=8), DIMENSION(mesh%dom_me*mesh%gauss%l_G,2), INTENT(IN), OPTIONAL :: opt_grad_div
204  REAL(KIND=8), DIMENSION(3*mesh%gauss%n_w) :: v145_loc, v236_loc
205  INTEGER, DIMENSION(3*mesh%gauss%n_w) :: idxm
206  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
207  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
208  REAL(KIND=8), DIMENSION(3,6) :: tensor_loc
209  REAL(KIND=8), DIMENSION(2) :: div_loc
210  REAL(KIND=8), DIMENSION(6) :: smb, f_tensor, f_div
211  INTEGER :: m, l, i, ni, index, nw, ix, ki, iglob, type, k
212  REAL(KIND=8) :: ray
213  vec :: vb_145, vb_236
214  petscerrorcode :: ierr
215 
216  CALL vecset(vb_145, 0.d0, ierr)
217  CALL vecset(vb_236, 0.d0, ierr)
218 
219  nw = mesh%gauss%n_w
220  index = 0
221  DO m = 1, mesh%dom_me
222  j_loc = mesh%jj(:,m)
223 
224  DO ki = 1, 3
225  DO ni = 1, nw
226  i = j_loc(ni)
227  iglob = vv_3_la%loc_to_glob(ki,i)
228  ix = (ki-1)*nw+ni
229  idxm(ix) = iglob-1
230  END DO
231  END DO
232 
233  v145_loc = 0.d0
234  v236_loc = 0.d0
235  DO l = 1, mesh%gauss%l_G
236  index = index + 1
237  dw_loc = mesh%gauss%dw(:,:,l,m)
238 
239  !===Compute radius of Gauss point
240  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
241 
242  IF (PRESENT(opt_tensor)) THEN
243  DO TYPE = 1, 6
244  DO k = 1, 3
245  tensor_loc(k,type) = opt_tensor(k,index,type)
246  END DO
247  END DO
248 
249  f_tensor(1) = (-mode*tensor_loc(1,4) + tensor_loc(2,3))/ray
250  f_tensor(2) = (mode*tensor_loc(1,3) + tensor_loc(2,4))/ray
251  f_tensor(3) = (-tensor_loc(1,3) - mode*tensor_loc(2,4))/ray
252  f_tensor(4) = (-tensor_loc(1,4) + mode*tensor_loc(2,3))/ray
253  f_tensor(5) = -tensor_loc(3,4)*mode/ray
254  f_tensor(6) = tensor_loc(3,3)*mode/ray
255  ELSE
256  f_tensor = 0.d0
257  END IF
258 
259  IF (PRESENT(opt_grad_div)) THEN
260  DO TYPE = 1, 2
261  div_loc(type) = opt_grad_div(index,type)
262  END DO
263 
264  f_div(1) = div_loc(1)/ray
265  f_div(2) = div_loc(2)/ray
266  f_div(3) = -mode*div_loc(2)/ray
267  f_div(4) = mode*div_loc(1)/ray
268  f_div(5) = 0.d0
269  f_div(6) = 0.d0
270  ELSE
271  f_div = 0.d0
272  END IF
273 
274  smb = (rhs_in(index,:)+f_tensor+f_div)*ray*mesh%gauss%rj(l,m)
275 
276  ki = 1
277  DO ni = 1, nw
278  ix = (ki-1)*nw+ni
279  v145_loc(ix) = v145_loc(ix) + mesh%gauss%ww(ni,l)*smb(1)
280  v236_loc(ix) = v236_loc(ix) + mesh%gauss%ww(ni,l)*smb(2)
281  END DO
282 
283  ki = 2
284  DO ni = 1, nw
285  ix = (ki-1)*nw+ni
286  v145_loc(ix) = v145_loc(ix) + mesh%gauss%ww(ni,l)*smb(4)
287  v236_loc(ix) = v236_loc(ix) + mesh%gauss%ww(ni,l)*smb(3)
288  END DO
289 
290  ki = 3
291  DO ni = 1, nw
292  ix = (ki-1)*nw+ni
293  v145_loc(ix) = v145_loc(ix) + mesh%gauss%ww(ni,l)*smb(5)
294  v236_loc(ix) = v236_loc(ix) + mesh%gauss%ww(ni,l)*smb(6)
295  END DO
296 
297  IF (PRESENT(opt_tensor)) THEN
298  ki = 1
299  DO ni = 1, nw
300  ix = (ki-1)*nw+ni
301  v145_loc(ix) = v145_loc(ix) &
302  + (dw_loc(1,ni)*tensor_loc(1,1) + dw_loc(2,ni)*tensor_loc(1,5))*ray*mesh%gauss%rj(l,m)
303  v236_loc(ix) = v236_loc(ix) &
304  + (dw_loc(1,ni)*tensor_loc(1,2) + dw_loc(2,ni)*tensor_loc(1,6))*ray*mesh%gauss%rj(l,m)
305  END DO
306  ki = 2
307  DO ni = 1, nw
308  ix = (ki-1)*nw+ni
309  v145_loc(ix) = v145_loc(ix) &
310  + (dw_loc(1,ni)*tensor_loc(2,2) + dw_loc(2,ni)*tensor_loc(2,6))*ray*mesh%gauss%rj(l,m)
311  v236_loc(ix) = v236_loc(ix) &
312  + (dw_loc(1,ni)*tensor_loc(2,1) + dw_loc(2,ni)*tensor_loc(2,5))*ray*mesh%gauss%rj(l,m)
313  END DO
314  ki = 3
315  DO ni = 1, nw
316  ix = (ki-1)*nw+ni
317  v145_loc(ix) = v145_loc(ix) &
318  + (dw_loc(1,ni)*tensor_loc(3,1) + dw_loc(2,ni)*tensor_loc(3,5))*ray*mesh%gauss%rj(l,m)
319  v236_loc(ix) = v236_loc(ix) &
320  + (dw_loc(1,ni)*tensor_loc(3,2) + dw_loc(2,ni)*tensor_loc(3,6))*ray*mesh%gauss%rj(l,m)
321  END DO
322  END IF
323 
324  IF (PRESENT(opt_grad_div)) THEN
325  ki = 1
326  DO ni = 1, nw
327  ix = (ki-1)*nw+ni
328  v145_loc(ix) = v145_loc(ix) &
329  + (dw_loc(1,ni)*div_loc(1))*ray*mesh%gauss%rj(l,m)
330  v236_loc(ix) = v236_loc(ix) &
331  + (dw_loc(1,ni)*div_loc(2))*ray*mesh%gauss%rj(l,m)
332  END DO
333  ki = 2 !nothing to do
334  ki = 3
335  DO ni = 1, nw
336  ix = (ki-1)*nw+ni
337  v145_loc(ix) = v145_loc(ix) &
338  + (dw_loc(2,ni)*div_loc(1))*ray*mesh%gauss%rj(l,m)
339  v236_loc(ix) = v236_loc(ix) &
340  + (dw_loc(2,ni)*div_loc(2))*ray*mesh%gauss%rj(l,m)
341  END DO
342  END IF
343 
344  ENDDO !l_G
345 
346  CALL vecsetvalues(vb_145, 3*nw, idxm, v145_loc, add_values, ierr)
347  CALL vecsetvalues(vb_236, 3*nw, idxm, v236_loc, add_values, ierr)
348 
349  ENDDO !dom_me
350 
351 
352  CALL vecassemblybegin(vb_145,ierr)
353  CALL vecassemblyend(vb_145,ierr)
354  CALL vecassemblybegin(vb_236,ierr)
355  CALL vecassemblyend(vb_236,ierr)
356  END SUBROUTINE rhs_3x3_art_comp
357 
358 
359 END MODULE rhs_para_assembling
def_type_mesh::petsc_csr_la
Definition: def_type_mesh.f90:18
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
my_util
Definition: my_util.f90:1
rhs_para_assembling
Definition: assembling_rhs.f90:1
rhs_para_assembling::rhs_3x3_art_comp
subroutine rhs_3x3_art_comp(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_grad_div)
Definition: assembling_rhs.f90:193
rhs_para_assembling::rhs_3x3
subroutine rhs_3x3(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_tensor_scaln_bdy)
Definition: assembling_rhs.f90:6
def_type_mesh
Definition: def_type_mesh.f90:4