subroutine process_element_hydraulic_linear_1(J, R, element, pressure, temperature, ice, porosity, &
properties, controls, actual_order)
implicit none
! --- arguments ---
type(type_crs), intent(inout) :: J
real(real64), intent(inout) :: R(:)
class(abst_mesh), intent(in), pointer :: element
type(type_variable), intent(in) :: pressure
type(type_variable), intent(in) :: temperature
type(type_variable), intent(in) :: ice
type(type_variable), intent(in) :: porosity
type(type_properties_manager), intent(in) :: properties
type(type_controls), intent(in) :: controls
integer(int32), intent(in) :: actual_order
! --- Local variables ---
integer(int32) :: index, num_nodes, num_gauss, material_id, il, jl, iG, iO
real(real64) :: weight, detJ
real(real64) :: dNdx_i, dNdy_i, dNdx_j, dNdy_j
real(real64) :: val, zeta
real(real64) :: dt
! --- Workspace variables ---
type(type_dense) :: CH_e
type(type_dense) :: KH_e
type(type_dense) :: J_e
real(real64), allocatable :: R_e(:)
! --- Physical quantities at Gauss points ---
type(type_state), allocatable :: state(:)
real(real64), dimension(:), pointer :: p_weight => null()
type(type_dp_vector_3d), dimension(:), pointer :: p_gauss => null()
real(real64), allocatable :: kflh(:)
real(real64), allocatable :: dot_ice(:)
integer(int32), dimension(:), pointer :: p_conn => null()
!---------------------------------------------------------------------------------------------------------------------------
! STEP 0: Initialize and obtain sizes
!---------------------------------------------------------------------------------------------------------------------------
num_nodes = element%get_num_nodes() !&
num_gauss = element%get_num_gauss() !&
material_id = element%get_group() !&
if (.not. controls%is_target(calc_hydraulic, material_id)) return
call CH_e%initialize_local(num_nodes)
call KH_e%initialize_local(num_nodes)
call J_e%initialize_local(num_nodes)
call allocate_array(R_e, num_nodes)
allocate (state(num_gauss))
call allocate_array(kflh, num_gauss)
call allocate_array(dot_ice, num_nodes)
dt = controls%time%get_dt()
p_weight => element%get_weight()
p_gauss => element%get_gauss()
p_conn => element%get_connectivity()
!---------------------------------------------------------------------------------------------------------------------------
! STEP 1: Compute the physical quantities at all Gauss points
!---------------------------------------------------------------------------------------------------------------------------
do iG = 1, num_gauss
state(iG)%temperature = element%lerp(p_gauss(iG), temperature%pre) !&
state(iG)%pressure = element%lerp(p_gauss(iG), pressure%pre) !&
state(iG)%porosity = element%lerp(p_gauss(iG), porosity%pre) !&
end do
call properties%calc_hydraulic(material_id, state, kflh)
do il = 1, num_nodes
dot_ice(il) = ice%dif(p_conn(il)) / dt
end do
!---------------------------------------------------------------------------------------------------------------------------
! STEP 2: Compute the element matrices CH_e and KH_e
!---------------------------------------------------------------------------------------------------------------------------
do iG = 1, num_gauss
weight = p_weight(iG) !&
detJ = element%jacobian_det(p_gauss(iG)) !&
zeta = state(iG)%density_ice / state(iG)%density_water - 1.0d0 !&
do il = 1, num_nodes
dNdx_i = ( element%jacobian(2, 2, p_gauss(iG)) * element%dpsi(il, 1, p_gauss(iG)) - & !&
element%jacobian(2, 1, p_gauss(iG)) * element%dpsi(il, 2, p_gauss(iG))) / detJ !&
dNdy_i = (-element%jacobian(1, 2, p_gauss(iG)) * element%dpsi(il, 1, p_gauss(iG)) + & !&
element%jacobian(1, 1, p_gauss(iG)) * element%dpsi(il, 2, p_gauss(iG))) / detJ !&
do jl = 1, num_nodes
dNdx_j = ( element%jacobian(2, 2, p_gauss(iG)) * element%dpsi(jl, 1, p_gauss(iG)) - & !&
element%jacobian(2, 1, p_gauss(iG)) * element%dpsi(jl, 2, p_gauss(iG))) / detJ !&
dNdy_j = (-element%jacobian(1, 2, p_gauss(iG)) * element%dpsi(jl, 1, p_gauss(iG)) + & !&
element%jacobian(1, 1, p_gauss(iG)) * element%dpsi(jl, 2, p_gauss(iG))) / detJ !&
val = element%psi(il, p_gauss(iG)) * element%psi(jl, p_gauss(iG)) * zeta * weight * detJ
call CH_e%add(il, jl, val)
val = (dNdx_i * dNdx_j + dNdy_i * dNdy_j) * kflh(iG) * weight * detJ
call KH_e%add(il, jl, val)
end do
end do
end do
!---------------------------------------------------------------------------------------------------------------------------
! STEP 3: Build the final local matrix (J_e) and vector (R_e)
!---------------------------------------------------------------------------------------------------------------------------
call add(0.0d0, KH_e, CH_e, J_e)
call gemv(-1.0d0, CH_e, dot_ice, 0.0d0, R_e)
!---------------------------------------------------------------------------------------------------------------------------
! STEP 4: Assemble the global matrix and vector
!---------------------------------------------------------------------------------------------------------------------------
do il = 1, num_nodes
R(p_conn(il)) = R(p_conn(il)) + R_e(il)
do jl = 1, num_nodes
call J%add(p_conn(il), p_conn(jl), J_e%val(il, jl))
end do
end do
!---------------------------------------------------------------------------------------------------------------------------
! STEP 5: Finalization
!---------------------------------------------------------------------------------------------------------------------------
call deallocate_array(kflh)
call deallocate_array(R_e)
call J_e%destroy()
call KH_e%destroy()
call CH_e%destroy()
deallocate (state)
end subroutine process_element_hydraulic_linear_1