← Back to branch summary

~fluidity-core/fluidity/sea-ice-branch

« back to all changes in this revision

Viewing changes to assemble/Assemble_CMC.F90

  • Committer: Simon Mouradian
  • Date: 2012-10-19 10:35:59 UTC
  • mfrom: (3520.32.371 fluidity)
  • Revision ID: simon.mouradian06@imperial.ac.uk-20121019103559-y36qa47phc69q8sc
merge from trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
assemble/Assemble_CMC.F90
172
172
 
173
173
    subroutine assemble_scaled_pressure_mass_matrix(state,scaled_pressure_mass_matrix)
174
174
 
175
 
      ! This routine assembles the scaled_pressure_mass_matrix. It is scaled
176
 
      ! by the inverse of viscosity.
 
175
      ! This routine assembles the scaled_pressure_mass_matrix at the
 
176
      ! quadrature points. It is scaled by the inverse of viscosity.
177
177
 
178
178
      type(state_type), intent(in) :: state   
179
 
 
180
179
      ! Scaled pressure mass matrix - already allocated in Momentum_Eq:
181
180
      type(csr_matrix), intent(inout) :: scaled_pressure_mass_matrix
182
 
 
183
181
      ! Pressure field:
184
182
      type(scalar_field), pointer :: pressure
185
183
      ! Viscosity tensor:
186
184
      type(tensor_field), pointer :: viscosity     
187
185
      ! Viscosity component:
188
186
      type(scalar_field) :: viscosity_component
189
 
      ! Inverse of Viscosity component
190
 
      type(scalar_field) :: inverse_viscosity_component
191
187
      ! Positions:
192
188
      type(vector_field), pointer :: positions
193
 
      
 
189
 
 
190
      ! Relevant declerations for mass matrix calculation:
 
191
      integer :: ele
 
192
      real, dimension(:), allocatable :: detwei
 
193
      type(element_type), pointer :: p_shape
 
194
      real, dimension(:,:), allocatable :: mass_matrix
 
195
      real, dimension(:), allocatable :: mu_gi
 
196
 
194
197
      ewrite(1,*) 'Entering assemble_scaled_pressure_mass_matrix'    
195
198
 
196
199
      ! Positions:
205
208
      ! Extract first component of viscosity tensor from full tensor:
206
209
      viscosity_component = extract_scalar_field(viscosity,1,1)
207
210
 
208
 
      ! Allocate memory for inverse viscosity component scalar field:
209
 
      call allocate(inverse_viscosity_component, viscosity_component%mesh, name="inverse_viscosity_component")
210
 
 
211
 
      ! Invert viscosity:
212
 
      call invert(viscosity_component,inverse_viscosity_component)
213
 
 
214
 
      ! Compute pressure mass matrix, scaled by inverse viscosity - note that the
215
 
      ! inverse viscosity is supplied as density here:
216
 
      call compute_mass(positions,pressure%mesh,scaled_pressure_mass_matrix,density=inverse_viscosity_component)
 
211
      ! Initialise and assemble scaled pressure mass matrix:
 
212
      allocate(detwei(ele_ngi(pressure, 1)), &
 
213
               mass_matrix(ele_loc(pressure, 1), ele_loc(pressure, 1)), &
 
214
               mu_gi(ele_ngi(viscosity_component, 1)))
 
215
 
 
216
      call zero(scaled_pressure_mass_matrix)
 
217
      do ele = 1, ele_count(pressure)
 
218
        p_shape => ele_shape(pressure, ele)
 
219
        mu_gi = ele_val_at_quad(viscosity_component, ele)
 
220
        call transform_to_physical(positions, ele, detwei=detwei)
 
221
        mass_matrix = shape_shape(p_shape, p_shape, detwei/mu_gi)
 
222
        call addto(scaled_pressure_mass_matrix, ele_nodes(pressure, ele),&
 
223
             ele_nodes(pressure, ele), mass_matrix)
 
224
      end do
217
225
 
218
226
      ewrite_minmax(scaled_pressure_mass_matrix)
219
227
 
220
 
      ! Deallocate inverse viscosity component scalar field:
221
 
      call deallocate(inverse_viscosity_component)
 
228
      deallocate(detwei, mass_matrix, mu_gi)
222
229
 
223
230
    end subroutine assemble_scaled_pressure_mass_matrix
224
231