SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
assembling_rhs.f90
Go to the documentation of this file.
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  IMPLICIT NONE
9  TYPE(mesh_type) :: mesh
10  TYPE(petsc_csr_la) :: vv_3_la
11  INTEGER, INTENT(IN) :: mode
12  REAL(KIND=8), DIMENSION(mesh%dom_me*mesh%gauss%l_G,6), INTENT(IN) :: rhs_in
13  REAL(KIND=8), DIMENSION(3,mesh%dom_me*mesh%gauss%l_G,6), INTENT(IN), OPTIONAL :: opt_tensor
14  REAL(KIND=8), DIMENSION(:,:), INTENT(IN), OPTIONAL :: opt_tensor_scaln_bdy
15  REAL(KIND=8), DIMENSION(3*mesh%gauss%n_w) :: v145_loc, v236_loc
16  INTEGER, DIMENSION(3*mesh%gauss%n_w) :: idxm
17  INTEGER, DIMENSION(3*mesh%gauss%n_ws) :: idxms
18  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
19  INTEGER, DIMENSION(mesh%gauss%n_ws) :: js_loc
20  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
21  REAL(KIND=8), DIMENSION(3,6) :: tensor_loc
22  REAL(KIND=8), DIMENSION(6) :: smb, f_tensor
23  INTEGER :: m, l, i, ni, index, nw, ix, ki, iglob, type, k, ms, ls, indexs, nws
24  REAL(KIND=8) :: ray
25 #include "petsc/finclude/petsc.h"
26  vec :: vb_145, vb_236
27  petscerrorcode :: ierr
28 
29  CALL vecset(vb_145, 0.d0, ierr)
30  CALL vecset(vb_236, 0.d0, ierr)
31 
32  nw = mesh%gauss%n_w
33  index = 0
34  DO m = 1, mesh%dom_me
35  j_loc = mesh%jj(:,m)
36 
37  DO ki = 1, 3
38  DO ni = 1, nw
39  i = j_loc(ni)
40  iglob = vv_3_la%loc_to_glob(ki,i)
41  ix = (ki-1)*nw+ni
42  idxm(ix) = iglob-1
43  END DO
44  END DO
45 
46  v145_loc = 0.d0
47  v236_loc = 0.d0
48  DO l = 1, mesh%gauss%l_G
49  index = index + 1
50  dw_loc = mesh%gauss%dw(:,:,l,m)
51 
52  !===Compute radius of Gauss point
53  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
54 
55  IF (present(opt_tensor)) THEN
56  DO TYPE = 1, 6
57  DO k = 1, 3
58  tensor_loc(k,type) = opt_tensor(k,index,type)
59  END DO
60  END DO
61 
62  f_tensor(1) = (-mode*tensor_loc(1,4) + tensor_loc(2,3))/ray
63  f_tensor(2) = (mode*tensor_loc(1,3) + tensor_loc(2,4))/ray
64  f_tensor(3) = (-tensor_loc(1,3) - mode*tensor_loc(2,4))/ray
65  f_tensor(4) = (-tensor_loc(1,4) + mode*tensor_loc(2,3))/ray
66  f_tensor(5) = -tensor_loc(3,4)*mode/ray
67  f_tensor(6) = tensor_loc(3,3)*mode/ray
68  ELSE
69  f_tensor = 0.d0
70  END IF
71 
72  smb = (rhs_in(index,:)+f_tensor)*ray*mesh%gauss%rj(l,m)
73 
74  ki = 1
75  DO ni = 1, nw
76  ix = (ki-1)*nw+ni
77  v145_loc(ix) = v145_loc(ix) + mesh%gauss%ww(ni,l)*smb(1)
78  v236_loc(ix) = v236_loc(ix) + mesh%gauss%ww(ni,l)*smb(2)
79  END DO
80 
81  ki = 2
82  DO ni = 1, nw
83  ix = (ki-1)*nw+ni
84  v145_loc(ix) = v145_loc(ix) + mesh%gauss%ww(ni,l)*smb(4)
85  v236_loc(ix) = v236_loc(ix) + mesh%gauss%ww(ni,l)*smb(3)
86  END DO
87 
88  ki = 3
89  DO ni = 1, nw
90  ix = (ki-1)*nw+ni
91  v145_loc(ix) = v145_loc(ix) + mesh%gauss%ww(ni,l)*smb(5)
92  v236_loc(ix) = v236_loc(ix) + mesh%gauss%ww(ni,l)*smb(6)
93  END DO
94 
95  IF (present(opt_tensor)) THEN
96  ki = 1
97  DO ni = 1, nw
98  ix = (ki-1)*nw+ni
99  v145_loc(ix) = v145_loc(ix) &
100  + (dw_loc(1,ni)*tensor_loc(1,1) + dw_loc(2,ni)*tensor_loc(1,5))*ray*mesh%gauss%rj(l,m)
101  v236_loc(ix) = v236_loc(ix) &
102  + (dw_loc(1,ni)*tensor_loc(1,2) + dw_loc(2,ni)*tensor_loc(1,6))*ray*mesh%gauss%rj(l,m)
103  END DO
104  ki = 2
105  DO ni = 1, nw
106  ix = (ki-1)*nw+ni
107  v145_loc(ix) = v145_loc(ix) &
108  + (dw_loc(1,ni)*tensor_loc(2,2) + dw_loc(2,ni)*tensor_loc(2,6))*ray*mesh%gauss%rj(l,m)
109  v236_loc(ix) = v236_loc(ix) &
110  + (dw_loc(1,ni)*tensor_loc(2,1) + dw_loc(2,ni)*tensor_loc(2,5))*ray*mesh%gauss%rj(l,m)
111  END DO
112  ki = 3
113  DO ni = 1, nw
114  ix = (ki-1)*nw+ni
115  v145_loc(ix) = v145_loc(ix) &
116  + (dw_loc(1,ni)*tensor_loc(3,1) + dw_loc(2,ni)*tensor_loc(3,5))*ray*mesh%gauss%rj(l,m)
117  v236_loc(ix) = v236_loc(ix) &
118  + (dw_loc(1,ni)*tensor_loc(3,2) + dw_loc(2,ni)*tensor_loc(3,6))*ray*mesh%gauss%rj(l,m)
119  END DO
120  END IF
121  ENDDO !l_G
122 
123  CALL vecsetvalues(vb_145, 3*nw, idxm, v145_loc, add_values, ierr)
124  CALL vecsetvalues(vb_236, 3*nw, idxm, v236_loc, add_values, ierr)
125 
126  ENDDO !dom_me
127 
128  IF (present(opt_tensor_scaln_bdy)) THEN
129  nws = mesh%gauss%n_ws
130  indexs = 0
131  DO ms = 1, mesh%dom_mes
132  js_loc = mesh%jjs(:,ms)
133 
134  DO ki = 1, 3
135  DO ni = 1, nws
136  i = js_loc(ni)
137  iglob = vv_3_la%loc_to_glob(ki,i)
138  ix = (ki-1)*nws+ni
139  idxms(ix) = iglob-1
140  END DO
141  END DO
142 
143  v145_loc = 0.d0
144  v236_loc = 0.d0
145  DO ls = 1, mesh%gauss%l_Gs
146  indexs = indexs + 1
147 
148  !===Compute radius of Gauss point
149  ray = sum(mesh%rr(1,js_loc)*mesh%gauss%wws(:,ls))
150 
151  !===Don't integrate on r=0
152  IF (ray.LT.1.d-10) cycle
153 
154  smb(:) = opt_tensor_scaln_bdy(indexs,:)*ray*mesh%gauss%rjs(ls,ms)
155 
156  ki = 1
157  DO ni = 1, nws
158  ix = (ki-1)*nws+ni
159  v145_loc(ix) = v145_loc(ix) + mesh%gauss%wws(ni,ls)*smb(1)
160  v236_loc(ix) = v236_loc(ix) + mesh%gauss%wws(ni,ls)*smb(2)
161  END DO
162 
163  ki = 2
164  DO ni = 1, nws
165  ix = (ki-1)*nws+ni
166  v145_loc(ix) = v145_loc(ix) + mesh%gauss%wws(ni,ls)*smb(4)
167  v236_loc(ix) = v236_loc(ix) + mesh%gauss%wws(ni,ls)*smb(3)
168  END DO
169 
170  ki = 3
171  DO ni = 1, nws
172  ix = (ki-1)*nws+ni
173  v145_loc(ix) = v145_loc(ix) + mesh%gauss%wws(ni,ls)*smb(5)
174  v236_loc(ix) = v236_loc(ix) + mesh%gauss%wws(ni,ls)*smb(6)
175  END DO
176 
177 
178  END DO
179  CALL vecsetvalues(vb_145, 3*nws, idxms, v145_loc, add_values, ierr)
180  CALL vecsetvalues(vb_236, 3*nws, idxms, v236_loc, add_values, ierr)
181  END DO
182  END IF
183 
184  CALL vecassemblybegin(vb_145,ierr)
185  CALL vecassemblyend(vb_145,ierr)
186  CALL vecassemblybegin(vb_236,ierr)
187  CALL vecassemblyend(vb_236,ierr)
188  END SUBROUTINE rhs_3x3
189 
190 
191 
192 END MODULE rhs_para_assembling
subroutine rhs_3x3(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_tensor_scaln_bdy)