1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
|
! Copyright (C) 2006 Imperial College London and others.
!
! Please see the AUTHORS file in the main source directory for a full list
! of copyright holders.
!
! Prof. C Pain
! Applied Modelling and Computation Group
! Department of Earth Science and Engineering
! Imperial College London
!
! amcgsoftware@imperial.ac.uk
!
! This library is free software; you can redistribute it and/or
! modify it under the terms of the GNU Lesser General Public
! License as published by the Free Software Foundation,
! version 2.1 of the License.
!
! This library is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
! Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public
! License along with this library; if not, write to the Free Software
! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
! USA
#include "fdebug.h"
module gradient_matrix_cg
use quadrature
use fields
use field_derivatives
use state_module
use futils
use fetools
use spud
use boundary_conditions
use global_parameters, only: OPTION_PATH_LEN
use transform_elements
use fldebug
use field_options, only: complete_field_path
implicit none
private
public :: assemble_gradient_matrix_cg
contains
subroutine assemble_gradient_matrix_cg(C_m, state, c_rhs, &
test_mesh, field, option_path, &
grad_mass, div_mass)
! inputs/outputs
! bucket full of fields
type(state_type), intent(inout) :: state
! the pressure gradient and compressible gradient matrices
type(block_csr_matrix), intent(inout) :: C_m
type(vector_field), intent(inout), optional :: c_rhs
type(mesh_type), intent(in) :: test_mesh
type(scalar_field), intent(inout) :: field
character(len=*), intent(in), optional :: option_path
type(csr_matrix), intent(inout), optional :: grad_mass, div_mass
! local
integer, dimension(:), pointer :: test_nodes, field_nodes
integer, dimension(:), allocatable :: test_nodes_bdy, field_nodes_bdy
real, dimension(:,:,:), allocatable :: ele_mat, ele_mat_bdy
type(element_type), pointer :: field_shape, test_shape
real, dimension(:,:,:), allocatable :: dfield_t
real, dimension(:,:,:), allocatable :: dtest_t
real, dimension(:), allocatable :: detwei, detwei_bdy
real, dimension(:,:), allocatable :: normal_bdy
! loop integers
integer :: ele, sele, dim
! pointer to coordinates
type(vector_field), pointer :: x
character(len=OPTION_PATH_LEN) :: l_option_path
! integrate by parts
logical :: integrate_by_parts
integer, dimension(:), allocatable :: field_bc_type
type(scalar_field) :: field_bc
real, dimension(:,:), allocatable :: div_mass_mat, grad_mass_mat
integer :: stat
! =============================================================
! Subroutine to construct the matrix CT_m (a.k.a. C1/2/3T).
! =============================================================
ewrite(2,*) 'In assemble_divergence_matrix_cg'
if(present(option_path)) then
l_option_path = trim(option_path)
else
l_option_path = trim(field%option_path)
end if
x=>extract_vector_field(state, "Coordinate")
integrate_by_parts=have_option(trim(complete_field_path(l_option_path, stat=stat))//&
&"/integrate_gradient_by_parts")
! Clear memory of arrays being designed
call zero(C_m)
if(present(c_rhs)) call zero(c_rhs)
if(present(div_mass)) call zero(div_mass)
if(present(grad_mass)) call zero(grad_mass)
allocate(dfield_t(ele_loc(field, 1), ele_ngi(field, 1), mesh_dim(field)), &
dtest_t(ele_loc(test_mesh, 1), ele_ngi(test_mesh, 1), mesh_dim(field)), &
ele_mat(mesh_dim(test_mesh), ele_loc(test_mesh, 1), ele_loc(field, 1)), &
detwei(ele_ngi(field, 1)), &
grad_mass_mat(ele_loc(test_mesh, 1), ele_loc(test_mesh, 1)), &
div_mass_mat(ele_loc(field, 1), ele_loc(field, 1)))
do ele=1, element_count(test_mesh)
test_nodes=>ele_nodes(test_mesh, ele)
field_nodes=>ele_nodes(field, ele)
test_shape=>ele_shape(test_mesh, ele)
field_shape=>ele_shape(field, ele)
if(integrate_by_parts) then
! transform the pressure derivatives into physical space
! (and get detwei)
call transform_to_physical(X, ele, test_shape, dshape=dtest_t,&
& detwei=detwei)
ele_mat = -dshape_shape(dtest_t, field_shape, detwei)
else
! transform the velociy derivatives into physical space
! (and get detwei)
call transform_to_physical(X, ele, field_shape, dshape=dfield_t,&
& detwei=detwei)
ele_mat = shape_dshape(test_shape, dfield_t, detwei)
end if
do dim = 1, mesh_dim(test_mesh)
call addto(c_m, dim, 1, test_nodes, field_nodes, ele_mat(dim,:,:))
end do
if(present(div_mass)) then
div_mass_mat = shape_shape(field_shape, field_shape, detwei)
call addto(div_mass, field_nodes, field_nodes, div_mass_mat)
end if
if(present(grad_mass)) then
grad_mass_mat = shape_shape(test_shape, test_shape, detwei)
call addto(grad_mass, test_nodes, test_nodes, grad_mass_mat)
end if
end do
if(integrate_by_parts) then
allocate(detwei_bdy(face_ngi(field, 1)), &
normal_bdy(mesh_dim(field), face_ngi(field, 1)))
allocate(field_nodes_bdy(field%mesh%faces%shape%loc))
allocate(test_nodes_bdy(test_mesh%faces%shape%loc))
allocate(ele_mat_bdy(mesh_dim(field), face_loc(test_mesh, 1), face_loc(field, 1)))
assert(surface_element_count(test_mesh)==surface_element_count(field))
allocate(field_bc_type(surface_element_count(field)))
call get_entire_boundary_condition(field, (/"weakdirichlet"/), field_bc, field_bc_type)
do sele = 1, surface_element_count(test_mesh)
test_shape=>face_shape(test_mesh, sele)
field_shape=>face_shape(field, sele)
test_nodes_bdy=face_global_nodes(test_mesh, sele)
field_nodes_bdy=face_global_nodes(field, sele)
call transform_facet_to_physical(X, sele, &
& detwei_f=detwei_bdy,&
& normal=normal_bdy)
ele_mat_bdy = shape_shape_vector(test_shape, field_shape, detwei_bdy, normal_bdy)
do dim = 1, mesh_dim(field)
if((field_bc_type(sele)==1).and.present(c_rhs)) then
call addto(c_rhs, dim, test_nodes_bdy, &
-matmul(ele_mat_bdy(dim,:,:), &
ele_val(field_bc, sele)))
else
call addto(c_m, dim, 1, test_nodes_bdy, field_nodes_bdy, &
ele_mat_bdy(dim,:,:))
end if
end do
end do
call deallocate(field_bc)
deallocate(field_bc_type)
deallocate(detwei_bdy, normal_bdy)
deallocate(test_nodes_bdy, field_nodes_bdy)
end if
deallocate(detwei)
end subroutine assemble_gradient_matrix_cg
end module gradient_matrix_cg
|