5 SUBROUTINE rhs_3x3(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_tensor_scaln_bdy)
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
25 #include "petsc/finclude/petsc.h"
27 petscerrorcode :: ierr
29 CALL vecset(vb_145, 0.d0, ierr)
30 CALL vecset(vb_236, 0.d0, ierr)
40 iglob = vv_3_la%loc_to_glob(ki,i)
48 DO l = 1, mesh%gauss%l_G
50 dw_loc = mesh%gauss%dw(:,:,l,m)
53 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
55 IF (present(opt_tensor))
THEN
58 tensor_loc(k,type) = opt_tensor(k,index,type)
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
72 smb = (rhs_in(index,:)+f_tensor)*ray*mesh%gauss%rj(l,m)
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)
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)
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)
95 IF (present(opt_tensor))
THEN
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)
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)
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)
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)
128 IF (present(opt_tensor_scaln_bdy))
THEN
129 nws = mesh%gauss%n_ws
131 DO ms = 1, mesh%dom_mes
132 js_loc = mesh%jjs(:,ms)
137 iglob = vv_3_la%loc_to_glob(ki,i)
145 DO ls = 1, mesh%gauss%l_Gs
149 ray = sum(mesh%rr(1,js_loc)*mesh%gauss%wws(:,ls))
152 IF (ray.LT.1.d-10) cycle
154 smb(:) = opt_tensor_scaln_bdy(indexs,:)*ray*mesh%gauss%rjs(ls,ms)
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)
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)
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)
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)
184 CALL vecassemblybegin(vb_145,ierr)
185 CALL vecassemblyend(vb_145,ierr)
186 CALL vecassemblybegin(vb_236,ierr)
187 CALL vecassemblyend(vb_236,ierr)
subroutine rhs_3x3(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_tensor_scaln_bdy)