Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_spline_fem_utilities_3d_helper.F90
Go to the documentation of this file.
1 
7 
8 module sll_m_spline_fem_utilities_3d_helper
9  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
10 #include "sll_assert.h"
11 #include "sll_errors.h"
12 #include "sll_memory.h"
13 #include "sll_working_precision.h"
14 
15  use sll_m_matrix_csr, only: &
16  sll_t_matrix_csr
17 
18  implicit none
19 
20  public :: &
21  sll_s_spline_fem_sparsity_mass3d, &
22  sll_s_spline_fem_sparsity_mass3d_clamped, &
23  sll_s_spline_fem_sparsity_mixedmass3d, &
24  sll_s_spline_fem_sparsity_mixedmass3d_clamped, &
25  assemble_mass3d, &
26  assemble_mass3d_clamped, &
27  assemble_mass3d_clamped_boundary, &
28  assemble_mixedmass3d, &
29  assemble_mixedmass3d_clamped_boundary, &
30  assemble_mixedmass3d_clamped
31 
32  private
33  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34 
35 contains
36 
37  !---------------------------------------------------------------------------!
39  subroutine sll_s_spline_fem_sparsity_mass3d( deg, n_cells, spmat )
40  sll_int32, intent( in ) :: deg(3)
41  sll_int32, intent( in ) :: n_cells(3)
42  type(sll_t_matrix_csr), intent( out ) :: spmat
43 
44  !local variables
45  sll_int32 :: n_nnz
46  sll_int32 :: n_nnz_per_row
47  sll_int32 :: ind, shift,row,i,j,k,k2
48 
49  ! Compute number of non-zero elements
50  n_nnz_per_row = (2*deg(1)+1)*(2*deg(2)+1)*(2*deg(3)+1)
51  n_nnz = n_nnz_per_row*n_cells(1)*n_cells(2)*n_cells(3)
52 
53  ! Create the csr matrix
54  call spmat%create( n_rows=n_cells(1)*n_cells(2)*n_cells(3), n_cols=n_cells(1)*n_cells(2)*n_cells(3), n_nnz=n_nnz )
55 
56  ! Set the row and column vector for the elements (row vector compressed)
57  ! There are n_nnz_per_row elements in each of the rows
58  spmat%arr_ia(1) = 1
59  do row = 2, n_cells(1)*n_cells(2)*n_cells(3)+1
60  spmat%arr_ia(row) = spmat%arr_ia(row-1) + n_nnz_per_row
61  end do
62 
63  sll_assert( spmat%arr_ia(n_cells(1)*n_cells(2)*n_cells(3)+1) == n_nnz+1 )
64 
65 
66  ind = 1
67  !initialise mass matrix
68  !loop over number of rows
69  do k = 1, n_cells(3)
70  if(k <= deg(3))then
71  do j = 1, n_cells(2)
72  do i = 1, n_cells(1)
73  do k2 = 1, k+deg(3)
74  shift = (k2-1)*n_cells(1)*n_cells(2)
75  call loop2( i, j, shift, deg, n_cells, ind, spmat )
76  end do
77  do k2 = n_cells(3)+k-deg(3), n_cells(3)
78  shift = (k2-1)*n_cells(1)*n_cells(2)
79  call loop2( i, j, shift, deg, n_cells, ind, spmat )
80  end do
81  sll_assert( ind == (i+(j-1)*n_cells(1)+(k-1)*n_cells(1)*n_cells(2))*n_nnz_per_row+1 )
82  end do
83  end do
84  elseif( k > deg(3) .and. k<= n_cells(3)-deg(3))then
85  do j = 1, n_cells(2)
86  do i = 1, n_cells(1)
87  do k2 = k-deg(3), k+deg(3)
88  shift = (k2-1)*n_cells(1)*n_cells(2)
89  call loop2( i, j, shift, deg, n_cells, ind, spmat )
90  end do
91  sll_assert( ind == (i+(j-1)*n_cells(1)+(k-1)*n_cells(1)*n_cells(2))*n_nnz_per_row+1 )
92  end do
93  end do
94  else
95  do j =1, n_cells(2)
96  do i =1, n_cells(1)
97  do k2 = 1, k+deg(3)-n_cells(3)
98  shift = (k2-1)*n_cells(1)*n_cells(2)
99  call loop2( i, j, shift, deg, n_cells, ind, spmat )
100  end do
101  do k2 = k-deg(3), n_cells(3)
102  shift = (k2-1)*n_cells(1)*n_cells(2)
103  call loop2( i, j, shift, deg, n_cells, ind, spmat )
104  end do
105  sll_assert( ind == (i+(j-1)*n_cells(1)+(k-1)*n_cells(1)*n_cells(2))*n_nnz_per_row+1 )
106  end do
107  end do
108  end if
109  end do
110  spmat%arr_a = 0.0_f64
111 
112  sll_assert( maxval(spmat%arr_ja) == spmat%n_cols )
113 
114  sll_assert( ind == spmat%n_nnz+1 )
115 
116  end subroutine sll_s_spline_fem_sparsity_mass3d
117 
118  !initialise sparse matrix in second direction
119  subroutine loop2( i, j, shift, deg, n_cells, ind, spmat )
120  sll_int32, intent( in ) :: i
121  sll_int32, intent( in ) :: j
122  sll_int32, intent( in ) :: shift
123  sll_int32, intent( in ) :: deg(3)
124  sll_int32, intent( in ) :: n_cells(3)
125  sll_int32, intent( inout ) :: ind
126  type(sll_t_matrix_csr), intent( inout ) :: spmat
127  !local variables
128  sll_int32 :: j2, shift1
129 
130  !in the first deg(2) rows the entries are splitted in the second direction
131  if(j<=deg(2))then
132  do j2 = 1, j+deg(2)
133  shift1 = (j2-1)*n_cells(1)+shift
134  call loop1( i, shift1, deg(1), n_cells(1), ind, spmat )
135  end do
136  do j2 = n_cells(2)+j-deg(2), n_cells(2)
137  shift1 = (j2-1)*n_cells(1)+shift
138  call loop1( i, shift1, deg(1), n_cells(1), ind, spmat )
139  end do
140 
141  !in rows: deg(2)+1-n_cells(2)-deg(2) we have consecutive entries in the second direction
142  elseif( j > deg(2) .and. j<= n_cells(2)-deg(2))then
143  do j2 = j-deg(2), j+deg(2)
144  shift1 = (j2-1)*n_cells(1)+shift
145  call loop1( i, shift1, deg(1), n_cells(1), ind, spmat )
146  end do
147 
148  !in the last deg(2) rows the entries are splitted in the second direction
149  else
150  do j2 = 1, j+deg(2)-n_cells(2)
151  shift1 = (j2-1)*n_cells(1)+shift
152  call loop1( i, shift1, deg(1), n_cells(1), ind, spmat )
153  end do
154  do j2 = j-deg(2), n_cells(2)
155  shift1 = (j2-1)*n_cells(1)+shift
156  call loop1(i,shift1, deg(1),n_cells(1),ind,spmat)
157  end do
158  end if
159 
160  end subroutine loop2
161 
162  !initialise sparse matrix in first direction
163  subroutine loop1( i, shift, deg, n_cells, ind, spmat )
164  sll_int32, intent( in ) :: i
165  sll_int32, intent( in ) :: shift
166  sll_int32, intent( in ) :: deg
167  sll_int32, intent( in ) :: n_cells
168  sll_int32, intent( inout ) :: ind
169  type(sll_t_matrix_csr), intent( inout ) :: spmat
170  !local variables
171  sll_int32 :: i2
172 
173  !in the first deg(1) rows the entries are splitted in the first direction
174  if(i<= deg)then
175  do i2 = 1, i+deg
176  spmat%arr_ja(ind) = i2+shift
177  ind = ind+1
178  end do
179  do i2 = n_cells+i-deg, n_cells
180  spmat%arr_ja(ind) = i2+shift
181  ind = ind+1
182  end do
183 
184  !in rows: deg(1)+1-n_cells(1)-deg(1) we have consecutive entries in the first direction
185  elseif( i > deg .and. i<= n_cells-deg) then
186  do i2 = i-deg, i+deg
187  spmat%arr_ja(ind) = i2+shift
188  ind = ind+1
189  end do
190 
191  !in the last deg(1) rows the entries are splitted in the first direction
192  else
193  do i2 = 1, i+deg-n_cells
194  spmat%arr_ja(ind) = i2+shift
195  ind = ind+1
196  end do
197  do i2 = i-deg, n_cells
198  spmat%arr_ja(ind) = i2+shift
199  ind = ind+1
200  end do
201  end if
202  end subroutine loop1
203 
204  !---------------------------------------------------------------------------!
206  subroutine sll_s_spline_fem_sparsity_mass3d_clamped( deg, n_cells, spmat )
207  sll_int32, intent( in ) :: deg(3)
208  sll_int32, intent( in ) :: n_cells(3)
209  type(sll_t_matrix_csr), intent( out ) :: spmat
210 
211  !local variables
212  sll_int32 :: n_nnz
213  sll_int32 :: n_nnz_per_row
214  sll_int32 :: ind, shift,row,i,j,k,k2
215 
216  ! Compute number of non-zero elements
217  n_nnz_per_row = (2*deg(1)+1)*(2*deg(2)+1)*(2*deg(3)+1)
218  n_nnz = (n_nnz_per_row*(n_cells(1)-deg(1))+(3*deg(1)**2+deg(1))*(2*deg(2)+1)*(2*deg(3)+1))*n_cells(2)*n_cells(3)
219 
220  ! Create the csr matrix
221  call spmat%create( n_rows=(n_cells(1)+deg(1))*n_cells(2)*n_cells(3), n_cols=(n_cells(1)+deg(1))*n_cells(2)*n_cells(3), n_nnz=n_nnz )
222 
223  ! Set the row and column vector for the elements (row vector compressed)
224  ! There are n_nnz_per_row elements in each of the rows
225  spmat%arr_ia(1) = 1
226  row=1
227  do k = 1, n_cells(3)
228  do j = 1, n_cells(2)
229  do i=2, deg(1)+1
230  row=row+1
231  spmat%arr_ia(row) = spmat%arr_ia(row-1) + (deg(1) + i - 1)*(2*deg(2)+1)*(2*deg(3)+1)
232  end do
233  do i = deg(1)+2, n_cells(1)+1
234  row=row+1
235  spmat%arr_ia(row) = spmat%arr_ia(row-1) + n_nnz_per_row
236  end do
237  do i=n_cells(1)+2, n_cells(1)+deg(1)+1
238  row=row+1
239  spmat%arr_ia(row) = spmat%arr_ia(row-1) + (2*deg(1) + n_cells(1)+2 - i)*(2*deg(2)+1)*(2*deg(3)+1)
240  end do
241  end do
242  end do
243 
244  sll_assert( row == (n_cells(1)+deg(1))*n_cells(2)*n_cells(3)+1)
245  sll_assert( spmat%arr_ia(row) == n_nnz+1 )
246 
247 
248  ind = 1
249  !initialise mass matrix
250  !loop over number of rows
251  do k=1, deg(3)
252  do j=1,n_cells(2)
253  do i=1,n_cells(1)+deg(1)
254  do k2 = 1, k+deg(3)
255  shift=(k2-1)*(n_cells(1)+deg(1))*n_cells(2)
256  call loop2_clamped( i, j, shift, deg, n_cells, ind, spmat )
257  end do
258 
259  do k2 = n_cells(3)+k-deg(3), n_cells(3)
260  shift=(k2-1)*(n_cells(1)+deg(1))*n_cells(2)
261  call loop2_clamped( i, j, shift, deg, n_cells, ind, spmat )
262  end do
263  end do
264  end do
265  end do
266  do k = deg(3)+1, n_cells(3)-deg(3)
267  do j=1,n_cells(2)
268  do i=1,n_cells(1)+deg(1)
269  do k2 = k-deg(3), k+deg(3)
270  shift=(k2-1)*(n_cells(1)+deg(1))*n_cells(2)
271  call loop2_clamped( i, j, shift, deg, n_cells, ind, spmat )
272  end do
273  end do
274  end do
275  end do
276  do k = n_cells(3)-deg(3)+1, n_cells(3)
277  do j=1,n_cells(2)
278  do i=1,n_cells(1)+deg(1)
279  do k2 = 1, k+deg(3)-n_cells(3)
280  shift=(k2-1)*(n_cells(1)+deg(1))*n_cells(2)
281  call loop2_clamped( i, j, shift, deg, n_cells, ind, spmat )
282  end do
283  do k2 = k-deg(3), n_cells(3)
284  shift=(k2-1)*(n_cells(1)+deg(1))*n_cells(2)
285  call loop2_clamped( i, j, shift, deg, n_cells, ind, spmat )
286  end do
287  end do
288  end do
289  end do
290  spmat%arr_a = 0.0_f64
291 
292  sll_assert( maxval(spmat%arr_ja) == spmat%n_cols )
293  sll_assert( ind == spmat%n_nnz+1 )
294 
295  end subroutine sll_s_spline_fem_sparsity_mass3d_clamped
296 
297  !initialise sparse matrix in second direction
298  subroutine loop2_clamped( i, j, shift, deg, n_cells, ind, spmat )
299  sll_int32, intent( in ) :: i
300  sll_int32, intent( in ) :: j
301  sll_int32, intent( in ) :: shift
302  sll_int32, intent( in ) :: deg(3)
303  sll_int32, intent( in ) :: n_cells(3)
304  sll_int32, intent( inout ) :: ind
305  type(sll_t_matrix_csr), intent( inout ) :: spmat
306  !local variables
307  sll_int32 :: j2, shift1
308 
309  !in the first deg(2) rows the entries are splitted in the second direction
310  if(j<=deg(2))then
311  do j2 = 1, j+deg(2)
312  shift1=(j2-1)*(n_cells(1)+deg(1))+shift
313  call loop1_clamped( i, shift1, deg(1), n_cells(1), ind, spmat )
314  end do
315  do j2 = n_cells(2)+j-deg(2), n_cells(2)
316  shift1=(j2-1)*(n_cells(1)+deg(1))+shift
317  call loop1_clamped( i, shift1, deg(1), n_cells(1), ind, spmat )
318  end do
319 
320  !in rows: deg(2)+1-n_cells(2)-deg(2) we have consecutive entries in the second direction
321  elseif( j > deg(2) .and. j<= n_cells(2)-deg(2) )then
322  do j2 = j-deg(2), j+deg(2)
323  shift1=(j2-1)*(n_cells(1)+deg(1))+shift
324  call loop1_clamped( i, shift1, deg(1), n_cells(1), ind, spmat )
325  end do
326 
327  !in the last deg(2) rows the entries are splitted in the second direction
328  elseif( j > n_cells(2)-deg(2) .and. j<= n_cells(2) )then
329  do j2 = 1, j+deg(2)-n_cells(2)
330  shift1=(j2-1)*(n_cells(1)+deg(1))+shift
331  call loop1_clamped( i, shift1, deg(1), n_cells(1), ind, spmat )
332  end do
333  do j2 = j-deg(2), n_cells(2)
334  shift1=(j2-1)*(n_cells(1)+deg(1))+shift
335  call loop1_clamped(i,shift1, deg(1),n_cells(1),ind,spmat)
336  end do
337  else
338  print*, 'error in loop2_clamped'
339  stop
340  end if
341 
342  end subroutine loop2_clamped
343 
344  !initialise sparse matrix in first direction
345  subroutine loop1_clamped( i, shift, deg, n_cells, ind, spmat )
346  sll_int32, intent( in ) :: i
347  sll_int32, intent( in ) :: shift
348  sll_int32, intent( in ) :: deg
349  sll_int32, intent( in ) :: n_cells
350  sll_int32, intent( inout ) :: ind
351  type(sll_t_matrix_csr), intent( inout ) :: spmat
352  !local variables
353  sll_int32 :: i2
354 
355  !in the first deg(1) rows we have between deg(1)+1 and 2*deg(1) entries
356  if(i <= deg)then
357  do i2 = 1, i+deg
358  spmat%arr_ja(ind) = i2+shift
359  ind = ind+1
360  end do
361  !in rows: deg(1)+1:n_cells(1) we have 2*deg(1)+1 consecutive entries in the first direction
362  elseif( i > deg .and. i <= n_cells) then
363  do i2 = i-deg, i+deg
364  spmat%arr_ja(ind) = i2+shift
365  ind = ind+1
366  end do
367  !in the last deg(1) rows we have between 2*deg(1) and deg(1)+1 entries
368  elseif( i > n_cells .and. i <= n_cells + deg ) then
369  do i2 = i-deg, n_cells+deg
370  spmat%arr_ja(ind) = i2+shift
371  ind = ind+1
372  end do
373  else
374  print*, 'error in loop1_clamped'
375  stop
376  end if
377 
378  end subroutine loop1_clamped
379 
380  !---------------------------------------------------------------------------!
382  subroutine sll_s_spline_fem_sparsity_mixedmass3d( deg1, deg2, n_cells, spmat )
383  sll_int32, intent( in ) :: deg1(3), deg2(3)
384  sll_int32, intent( in ) :: n_cells(3)
385  type(sll_t_matrix_csr), intent( out ) :: spmat
386 
387  !local variables
388  sll_int32 :: n_nnz
389  sll_int32 :: n_nnz_per_row
390  sll_int32 :: ind, shift,row,i,j,k,l
391 
392  ! Compute number of non-zero elements
393  n_nnz_per_row = (deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)*(deg1(3)+deg2(3)+1)
394  n_nnz = n_nnz_per_row*n_cells(1)*n_cells(2)*n_cells(3)
395  ! Create the csr matrix
396  call spmat%create( n_rows=n_cells(1)*n_cells(2)*n_cells(3), n_cols=n_cells(1)*n_cells(2)*n_cells(3), n_nnz=n_nnz )
397 
398  ! Set the row and column vector for the elements (row vector compressed)
399  ! There are n_nnz_per_row elements in each of the rows
400  spmat%arr_ia(1) = 1
401  do row = 2, n_cells(1)*n_cells(2)*n_cells(3)+1
402  spmat%arr_ia(row) = spmat%arr_ia(row-1) + n_nnz_per_row
403  end do
404 
405  sll_assert( spmat%arr_ia(n_cells(1)*n_cells(2)*n_cells(3)+1) == n_nnz+1 )
406 
407 
408  ind = 1
409  !initialise mass matrix
410  !loop over number of rows
411  do k = 1 ,n_cells(3)
412  if(k <= deg2(3)) then
413  do j = 1, n_cells(2)
414  do i = 1, n_cells(1)
415  do l = 1, k+deg1(3)
416  shift = (l-1)*n_cells(1)*n_cells(2)
417  call mloop2( i, j, shift, deg1, deg2, n_cells, ind, spmat )
418  end do
419  do l = n_cells(3)+k-deg2(3), n_cells(3)
420  shift = (l-1)*n_cells(1)*n_cells(2)
421  call mloop2( i, j, shift, deg1, deg2, n_cells, ind, spmat )
422  end do
423  sll_assert( ind == (i+(j-1)*n_cells(1)+(k-1)*n_cells(1)*n_cells(2))*n_nnz_per_row+1 )
424  end do
425  end do
426  elseif( k > deg2(3) .and. k <= n_cells(3)-deg1(3)) then
427  do j = 1, n_cells(2)
428  do i = 1, n_cells(1)
429  do l = k-deg2(3), k+deg1(3)
430  shift = (l-1)*n_cells(1)*n_cells(2)
431  call mloop2( i, j, shift, deg1, deg2, n_cells, ind, spmat )
432  end do
433  sll_assert( ind == (i+(j-1)*n_cells(1)+(k-1)*n_cells(1)*n_cells(2))*n_nnz_per_row+1 )
434  end do
435  end do
436  else
437  do j = 1, n_cells(2)
438  do i = 1, n_cells(1)
439  do l = 1, k+deg1(3)-n_cells(3)
440  shift = (l-1)*n_cells(1)*n_cells(2)
441  call mloop2( i, j, shift, deg1, deg2, n_cells, ind, spmat )
442  end do
443  do l = k-deg2(3), n_cells(3)
444  shift = (l-1)*n_cells(1)*n_cells(2)
445  call mloop2( i, j, shift, deg1, deg2, n_cells, ind, spmat )
446  end do
447  sll_assert( ind == (i+(j-1)*n_cells(1)+(k-1)*n_cells(1)*n_cells(2))*n_nnz_per_row+1 )
448  end do
449  end do
450  end if
451  end do
452  spmat%arr_a = 0.0_f64
453 
454  sll_assert( maxval(spmat%arr_ja) == spmat%n_cols )
455  sll_assert( ind == spmat%n_nnz+1 )
456 
457  end subroutine sll_s_spline_fem_sparsity_mixedmass3d
458 
459  !initialise sparse matrix in second direction
460  subroutine mloop2( i, j, shift, deg1, deg2, n_cells, ind, spmat )
461  sll_int32, intent( in ) :: i
462  sll_int32, intent( in ) :: j
463  sll_int32, intent( in ) :: shift
464  sll_int32, intent( in ) :: deg1(3), deg2(3)
465  sll_int32, intent( in ) :: n_cells(3)
466  sll_int32, intent( inout ) :: ind
467  type(sll_t_matrix_csr), intent( inout ) :: spmat
468  !local variables
469  sll_int32 :: l, shift1
470 
471  !in the first deg2(2) rows the entries are splitted in the second direction
472  if(j <= deg2(2))then
473  do l = 1, j+deg1(2)
474  shift1 = (l-1)*n_cells(1)+shift
475  call mloop1( i, shift1, deg1(1), deg2(1), n_cells(1), ind, spmat )
476  end do
477  do l = n_cells(2)+j-deg2(2), n_cells(2)
478  shift1 = (l-1)*n_cells(1)+shift
479  call mloop1( i, shift1, deg1(1), deg2(1), n_cells(1), ind, spmat )
480  end do
481 
482  !in rows: deg2(2)+1-n_cells(2)-deg1(2) we have consecutive entries in the second direction
483  elseif( j > deg2(2) .and. j <= n_cells(2)-deg1(2))then
484  do l = j-deg2(2), j+deg1(2)
485  shift1 = (l-1)*n_cells(1)+shift
486  call mloop1( i, shift1, deg1(1), deg2(1), n_cells(1), ind, spmat )
487  end do
488 
489  !in the last deg1(2) rows the entries are splitted in the second direction
490  else
491  do l = 1, j+deg1(2)-n_cells(2)
492  shift1 = (l-1)*n_cells(1)+shift
493  call mloop1( i, shift1, deg1(1), deg2(1), n_cells(1), ind, spmat )
494  end do
495  do l = j-deg2(2), n_cells(2)
496  shift1 = (l-1)*n_cells(1)+shift
497  call mloop1(i, shift1, deg1(1), deg2(1), n_cells(1), ind,spmat)
498  end do
499  end if
500 
501  end subroutine mloop2
502 
503  !initialise sparse matrix in first direction
504  subroutine mloop1( i, shift, deg1, deg2, n_cells, ind, spmat )
505  sll_int32, intent( in ) :: i
506  sll_int32, intent( in ) :: shift
507  sll_int32, intent( in ) :: deg1, deg2
508  sll_int32, intent( in ) :: n_cells
509  sll_int32, intent( inout ) :: ind
510  type(sll_t_matrix_csr), intent( inout ) :: spmat
511  !local variables
512  sll_int32 :: l
513 
514  !in the first deg2(1) rows the entries are splitted in the first direction
515  if(i <= deg2)then
516  do l = 1, i+deg1
517  spmat%arr_ja(ind) = l+shift
518  ind = ind+1
519  end do
520  do l = n_cells+i-deg2, n_cells
521  spmat%arr_ja(ind) = l+shift
522  ind = ind+1
523  end do
524 
525  !in rows: deg2(1)+1 to n_cells(1)-deg1(1) we have consecutive entries in the first direction
526  elseif( i > deg2 .and. i<= n_cells-deg1) then
527  do l = i-deg2, i+deg1
528  spmat%arr_ja(ind) = l+shift
529  ind = ind+1
530  end do
531 
532  !in the last deg1(1) rows the entries are splitted in the first direction
533  elseif(i > n_cells-deg1 .and. i <= n_cells) then
534  do l = 1, i+deg1-n_cells
535  spmat%arr_ja(ind) = l+shift
536  ind = ind+1
537  end do
538  do l = i-deg2, n_cells
539  spmat%arr_ja(ind) = l+shift
540  ind = ind+1
541  end do
542  end if
543  end subroutine mloop1
544 
545  !---------------------------------------------------------------------------!
547  subroutine sll_s_spline_fem_sparsity_mixedmass3d_clamped( deg1, deg2, n_cells, spmat )
548  sll_int32, intent( in ) :: deg1(3), deg2(3)
549  sll_int32, intent( in ) :: n_cells(3)
550  type(sll_t_matrix_csr), intent( out ) :: spmat
551  !local variables
552  sll_int32 :: n_nnz
553  sll_int32 :: n_nnz_per_row
554  sll_int32 :: ind, shift,row,i,j,k,l
555 
556  ! Compute number of non-zero elements
557  n_nnz_per_row = (deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)*(deg1(3)+deg2(3)+1)
558  n_nnz = (n_nnz_per_row*(n_cells(1)-deg1(1))+deg1(1)*(2*deg2(1)+deg1(1)+1)*(deg1(2)+deg2(2)+1)*(deg1(3)+deg2(3)+1))*n_cells(2)*n_cells(3)
559  ! Create the csr matrix
560  call spmat%create( n_rows=(n_cells(1)+deg1(1))*n_cells(2)*n_cells(3), n_cols=(n_cells(1)+deg2(1))*n_cells(2)*n_cells(3), n_nnz=n_nnz )
561 
562  ! Set the row and column vector for the elements (row vector compressed)
563  ! There are n_nnz_per_row elements in each of the rows
564  spmat%arr_ia(1) = 1
565  row=1
566  do k = 1, n_cells(3)
567  do j = 1, n_cells(2)
568  do i=2, deg1(1)+1
569  row=row+1
570  spmat%arr_ia(row) = spmat%arr_ia(row-1) + (deg2(1) + i - 1)*(deg1(2)+deg2(2)+1)*(deg1(3)+deg2(3)+1)
571  end do
572  do i = deg1(1)+2, n_cells(1)+1
573  row=row+1
574  spmat%arr_ia(row) = spmat%arr_ia(row-1) + n_nnz_per_row
575  end do
576  do i=n_cells(1)+2, n_cells(1)+deg1(1)+1
577  row=row+1
578  spmat%arr_ia(row) = spmat%arr_ia(row-1) + (deg1(1)+deg2(1) + n_cells(1)+2 - i)*(deg1(2)+deg2(2)+1)*(deg1(3)+deg2(3)+1)
579  end do
580  end do
581  end do
582 
583  sll_assert( row == (n_cells(1)+deg1(1))*n_cells(2)*n_cells(3)+1)
584  sll_assert( spmat%arr_ia(row) == n_nnz+1 )
585 
586 
587  ind = 1
588  !initialise mass matrix
589  !loop over number of rows
590  do k=1,n_cells(3)
591  if(k<=deg2(3)) then
592  do j=1,n_cells(2)
593  do i=1,n_cells(1)+deg1(1)
594  do l = 1, k+deg1(3)
595  shift=(l-1)*(n_cells(1)+deg2(1))*n_cells(2)
596  call mloop2_clamped( i, j, shift, deg1, deg2, n_cells, ind, spmat )
597  end do
598  do l = n_cells(3)+k-deg2(3), n_cells(3)
599  shift=(l-1)*(n_cells(1)+deg2(1))*n_cells(2)
600  call mloop2_clamped( i, j, shift, deg1, deg2, n_cells, ind, spmat )
601  end do
602  end do
603  end do
604  elseif( k > deg2(3) .and. k<= n_cells(3)-deg1(3)) then
605  do j=1,n_cells(2)
606  do i=1,n_cells(1)+deg1(1)
607  do l = k-deg2(3), k+deg1(3)
608  shift=(l-1)*(n_cells(1)+deg2(1))*n_cells(2)
609  call mloop2_clamped( i, j, shift, deg1, deg2, n_cells, ind, spmat )
610  end do
611  end do
612  end do
613  else
614  do j=1,n_cells(2)
615  do i=1, n_cells(1)+deg1(1)
616  do l = 1, k+deg1(3)-n_cells(3)
617  shift=(l-1)*(n_cells(1)+deg2(1))*n_cells(2)
618  call mloop2_clamped( i, j, shift, deg1, deg2, n_cells, ind, spmat )
619  end do
620  do l = k-deg2(3), n_cells(3)
621  shift=(l-1)*(n_cells(1)+deg2(1))*n_cells(2)
622  call mloop2_clamped( i, j, shift, deg1, deg2, n_cells, ind, spmat )
623  end do
624  end do
625  end do
626  end if
627  end do
628  spmat%arr_a = 0.0_f64
629 
630  sll_assert( maxval(spmat%arr_ja) == spmat%n_cols )
631  sll_assert( ind == spmat%n_nnz+1 )
632 
633  end subroutine sll_s_spline_fem_sparsity_mixedmass3d_clamped
634 
635  !initialise sparse mass matrix in second direction
636  subroutine mloop2_clamped( i, j, shift, deg1, deg2, n_cells, ind, spmat )
637  sll_int32, intent( in ) :: i
638  sll_int32, intent( in ) :: j
639  sll_int32, intent( in ) :: shift
640  sll_int32, intent( in ) :: deg1(3), deg2(3)
641  sll_int32, intent( in ) :: n_cells(3)
642  sll_int32, intent( inout ) :: ind
643  type(sll_t_matrix_csr), intent( inout ) :: spmat
644  !local variables
645  sll_int32 :: l, shift1
646 
647  !in the first deg2(2) rows the entries are splitted in the second direction
648  if(j<=deg2(2))then
649  do l = 1, j+deg1(2)
650  shift1=(l-1)*(n_cells(1)+deg2(1))+shift
651  call mloop1_clamped( i, shift1, deg1(1), deg2(1), n_cells(1), ind, spmat )
652  end do
653  do l = n_cells(2)+j-deg2(2), n_cells(2)
654  shift1=(l-1)*(n_cells(1)+deg2(1))+shift
655  call mloop1_clamped( i, shift1, deg1(1), deg2(1), n_cells(1), ind, spmat )
656  end do
657 
658  !in rows: deg2(2)+1-n_cells(2)-deg1(2) we have consecutive entries in the second direction
659  elseif( j > deg2(2) .and. j<= n_cells(2)-deg1(2))then
660  do l = j-deg2(2), j+deg1(2)
661  shift1=(l-1)*(n_cells(1)+deg2(1))+shift
662  call mloop1_clamped( i, shift1, deg1(1), deg2(1), n_cells(1), ind, spmat )
663  end do
664 
665  !in the last deg1(2) rows the entries are splitted in the second direction
666  else
667  do l = 1, j+deg1(2)-n_cells(2)
668  shift1=(l-1)*(n_cells(1)+deg2(1))+shift
669  call mloop1_clamped( i, shift1, deg1(1), deg2(1), n_cells(1), ind, spmat )
670  end do
671  do l = j-deg2(2), n_cells(2)
672  shift1=(l-1)*(n_cells(1)+deg2(1))+shift
673  call mloop1_clamped(i,shift1, deg1(1), deg2(1),n_cells(1),ind,spmat)
674  end do
675  end if
676 
677  end subroutine mloop2_clamped
678 
679  !initialise sparse mass matrix in first direction
680  subroutine mloop1_clamped( i, shift, deg1, deg2, n_cells, ind, spmat )
681  sll_int32, intent( in ) :: i
682  sll_int32, intent( in ) :: shift
683  sll_int32, intent( in ) :: deg1, deg2
684  sll_int32, intent( in ) :: n_cells
685  sll_int32, intent( inout ) :: ind
686  type(sll_t_matrix_csr), intent( inout ) :: spmat
687  !local variables
688  sll_int32 :: l
689 
690  !in the first deg1(1) rows we have between deg2+1 and deg2+deg1 entries
691  if(i<= deg1)then
692  do l = 1, i+deg2
693  spmat%arr_ja(ind) = l+shift
694  ind = ind+1
695  end do
696  !in rows: deg1+1 to n_cells we have deg1+deg2+1 consecutive entries in the first direction
697  elseif( i > deg1 .and. i<= n_cells) then
698  do l = i-deg1, i+deg2
699  spmat%arr_ja(ind) = l+shift
700  ind = ind+1
701  end do
702  !in the last deg1 rows we have between deg2+deg1 and deg2+1 entries
703  elseif( i > n_cells .and. i <= n_cells + deg1 ) then
704  do l = i-deg1, n_cells+deg2
705  spmat%arr_ja(ind) = l+shift
706  ind = ind+1
707  end do
708  end if
709  end subroutine mloop1_clamped
710 
711 
712  !---------------------------------------------------------------------------!
714  subroutine assemble_mass3d ( deg, n_cells, mass_line, matrix, row, ind )
715  sll_int32, intent( in ) :: deg(3)
716  sll_int32, intent( in ) :: n_cells(3)
717  sll_real64, intent( in ) :: mass_line((2*deg(1)+1)*(2*deg(2)+1)*(2*deg(3)+1))
718  type(sll_t_matrix_csr), intent( inout ) :: matrix
719  sll_int32, intent( in ) :: row(3)
720  sll_int32, intent( inout ) :: ind
721  !local variables
722  sll_int32 :: column, shift
723 
724  !assemble massmatrix in third direction
725  ! For the first deg(3) rows we need to put the first part to the back due to periodic boundaries
726  if(row(3) <= deg(3)) then
727  do column = 2-row(3)+deg(3), deg(3)
728  shift = (column-1)*(2*deg(1)+1)*(2*deg(2)+1)
729  call assemble2( deg, n_cells, mass_line, matrix, row, ind, shift )
730  end do
731  do column = deg(3)+1,2*deg(3)+1
732  shift = (column-1)*(2*deg(1)+1)*(2*deg(2)+1)
733  call assemble2( deg, n_cells, mass_line, matrix, row, ind, shift )
734  end do
735 
736  do column = 1, deg(3)-row(3)+1
737  shift=(column-1)*(2*deg(1)+1)*(2*deg(2)+1)
738  call assemble2( deg, n_cells, mass_line, matrix, row, ind, shift )
739  end do
740  else if(row(3) >= deg(3)+1 .and. row(3)<= n_cells(3)-deg(3)) then
741  do column = 1, 2*deg(3)+1
742  shift = (column-1)*(2*deg(1)+1)*(2*deg(2)+1)
743  call assemble2( deg, n_cells, mass_line, matrix, row, ind, shift )
744  end do
745 
746  ! For the last deg(3) rows, we need to put the second part to the front due to periodic bounaries
747  else if(row(3) >= n_cells(3)-deg(3)+1 .and. row(3) <= n_cells(3)) then
748 
749  do column = n_cells(3)-row(3)+deg(3)+2, 2*deg(3)+1
750  shift = (column-1)*(2*deg(1)+1)*(2*deg(2)+1)
751  call assemble2( deg, n_cells, mass_line, matrix, row, ind, shift )
752  end do
753 
754  do column = 1, deg(3)+1
755  shift = (column-1)*(2*deg(1)+1)*(2*deg(2)+1)
756  call assemble2( deg, n_cells, mass_line, matrix, row, ind, shift )
757  end do
758 
759  do column = 2+deg(3), n_cells(3)-row(3)+deg(3)+1
760  shift = (column-1)*(2*deg(1)+1)*(2*deg(2)+1)
761  call assemble2( deg, n_cells, mass_line, matrix, row, ind, shift )
762  end do
763  else
764  sll_error('assemble','error in row in assemble_mass')
765  end if
766 
767  sll_assert( ind == (row(1)+(row(2)-1)*n_cells(1)+(row(3)-1)*n_cells(1)*n_cells(2))*(2*deg(1)+1)*(2*deg(2)+1)*(2*deg(3)+1)+1)
768 
769  end subroutine assemble_mass3d
770 
771  !assemble mass matrix in second direction
772  subroutine assemble2( deg, n_cells, mass_line, matrix, row, ind, shift )
773  sll_int32, intent( in ) :: deg(3)
774  sll_int32, intent( in ) :: n_cells(3)
775  sll_real64, intent( in ) :: mass_line((2*deg(1)+1)*(2*deg(2)+1)*(2*deg(3)+1))
776  type(sll_t_matrix_csr), intent( inout ) :: matrix
777  sll_int32, intent( in ) :: row(3)
778  sll_int32, intent( inout ) :: ind
779  sll_int32, intent( in ) :: shift
780  !local variables
781  sll_int32 :: column, shift1
782  ! For the first deg(2) rows we need to put the first part to the back due to periodic boundaries
783  if(row(2) <= deg(2)) then
784  do column = 2-row(2)+deg(2), deg(2)
785  shift1 = (column-1)*(2*deg(1)+1)+shift
786  call assemble1( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
787  end do
788  do column = deg(2)+1, 2*deg(2)+1
789  shift1 = (column-1)*(2*deg(1)+1)+shift
790  call assemble1( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
791  end do
792 
793  do column = 1, deg(2)-row(2)+1
794  shift1 = (column-1)*(2*deg(1)+1)+shift
795  call assemble1( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
796  end do
797  else if(row(2) >= deg(2)+1 .and. row(2)<= n_cells(2)-deg(2)) then
798  do column = 1, 2*deg(2)+1
799  shift1 = (column-1)*(2*deg(1)+1)+shift
800  call assemble1( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
801  end do
802 
803  ! For the last deg(2) rows, we need to put the second part to the front due to periodic boundaries
804  else if(row(2) >= n_cells(2)-deg(2)+1 .and. row(2) <= n_cells(2)) then
805 
806  do column = n_cells(2)-row(2)+deg(2)+2, 2*deg(2)+1
807  shift1 = (column-1)*(2*deg(1)+1)+shift
808  call assemble1( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
809  end do
810 
811  do column = 1, deg(2)+1
812  shift1 = (column-1)*(2*deg(1)+1)+shift
813  call assemble1(deg, n_cells(1), mass_line, matrix, row(1), ind, shift1)
814  end do
815 
816  do column = 2+deg(2), n_cells(2)-row(2)+deg(2)+1
817  shift1 = (column-1)*(2*deg(1)+1)+shift
818  call assemble1( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
819  end do
820  else
821  sll_error('assemble','error in row in assemble_mass')
822  end if
823 
824  end subroutine assemble2
825 
826  !assemble mass matrix in first direction
827  subroutine assemble1( deg, n_cells, mass_line, matrix, row, ind, shift )
828  sll_int32, intent( in ) :: deg(3)
829  sll_int32, intent( in ) :: n_cells
830  sll_real64, intent( in ) :: mass_line((2*deg(1)+1)*(2*deg(2)+1)*(2*deg(3)+1))
831  type(sll_t_matrix_csr), intent( inout ) :: matrix
832  sll_int32, intent( in ) :: row
833  sll_int32, intent( inout ) :: ind
834  sll_int32, intent( in ) :: shift
835  !local variables
836  sll_int32 :: column
837 
838  ! For the first deg(1) rows we need to put the first part to the back due to periodic boundaries
839  if(row <= deg(1)) then
840  do column = 2-row+deg(1), deg(1)
841  matrix%arr_a(1,ind) = mass_line(column+shift)
842 
843  ind = ind+1
844  end do
845  matrix%arr_a(1,ind:ind+deg(1)) = mass_line(deg(1)+1+shift:2*deg(1)+1+shift)
846  ind = ind+deg(1)+1
847 
848  do column = 1, deg(1)-row+1
849  matrix%arr_a(1,ind) = mass_line(column+shift)
850  ind = ind+1
851  end do
852  else if(row > deg(1) .and. row <= n_cells-deg(1)) then
853  matrix%arr_a(1,ind:ind+2*deg(1)) = mass_line(1+shift:2*deg(1)+1+shift)
854 
855  ind = ind+2*deg(1)+1
856  ! For the last deg(1) rows, we need to put the second part to the front due to periodic boundaries
857  else if(row >= n_cells-deg(1)+1 .and. row <= n_cells) then
858 
859  do column = n_cells-row+deg(1)+2, 2*deg(1)+1
860  matrix%arr_a(1,ind) = mass_line(column+shift)
861  ind = ind+1
862  end do
863 
864  matrix%arr_a(1,ind:ind+deg(1)) = mass_line(1+shift:deg(1)+1+shift)
865  ind = ind+deg(1)+1
866 
867  do column = 2+deg(1), n_cells-row+deg(1)+1
868  matrix%arr_a(1,ind) = mass_line(column+shift)
869  ind = ind+1
870  end do
871  else
872  sll_error('assemble','error in row in assemble_mass')
873  end if
874 
875  end subroutine assemble1
876 
877  !---------------------------------------------------------------------------!
879  subroutine assemble_mass3d_clamped ( deg, n_cells, mass_line, matrix, row, ind )
880  sll_int32, intent( in ) :: deg(3)
881  sll_int32, intent( in ) :: n_cells(3)
882  sll_real64, intent( in ) :: mass_line((2*deg(1)+1)*(2*deg(2)+1)*(2*deg(3)+1))
883  type(sll_t_matrix_csr), intent( inout ) :: matrix
884  sll_int32, intent( in ) :: row(3)
885  sll_int32, intent( inout ) :: ind
886  !local variables
887  sll_int32 :: column, shift, begin
888 
889  begin= ind
890 
891  !assemble massmatrix in third direction
892  ! For the first deg(3) rows we need to put the first part to the back due to periodic boundaries
893  if(row(3)<=deg(3)) then
894  do column = 2-row(3)+deg(3), deg(3)
895  shift=(column-1)*(2*deg(1)+1)*(2*deg(2)+1)
896  call assemble2_clamped( deg, n_cells, mass_line, matrix, row, ind, shift )
897  end do
898  do column= deg(3)+1,2*deg(3)+1
899  shift=(column-1)*(2*deg(1)+1)*(2*deg(2)+1)
900  call assemble2_clamped( deg, n_cells, mass_line, matrix, row, ind, shift )
901  end do
902 
903  do column = 1, deg(3)-row(3)+1
904  shift=(column-1)*(2*deg(1)+1)*(2*deg(2)+1)
905  call assemble2_clamped( deg, n_cells, mass_line, matrix, row, ind, shift )
906  end do
907  else if(row(3) >= deg(3)+1 .and. row(3)<= n_cells(3)-deg(3)) then
908  do column=1, 2*deg(3)+1
909  shift=(column-1)*(2*deg(1)+1)*(2*deg(2)+1)
910  call assemble2_clamped( deg, n_cells, mass_line, matrix, row, ind, shift )
911  end do
912 
913  ! For the last deg(3) rows, we need to put the second part to the front due to periodic boundaries
914  else if(row(3) >= n_cells(3)-deg(3)+1 .and. row(3) <= n_cells(3)) then
915 
916  do column = n_cells(3)-row(3)+deg(3)+2, 2*deg(3)+1
917  shift=(column-1)*(2*deg(1)+1)*(2*deg(2)+1)
918  call assemble2_clamped( deg, n_cells, mass_line, matrix, row, ind, shift )
919  end do
920 
921  do column=1, deg(3)+1
922  shift=(column-1)*(2*deg(1)+1)*(2*deg(2)+1)
923  call assemble2_clamped( deg, n_cells, mass_line, matrix, row, ind, shift )
924  end do
925 
926  do column = 2+deg(3), n_cells(3)-row(3)+deg(3)+1
927  shift=(column-1)*(2*deg(1)+1)*(2*deg(2)+1)
928  call assemble2_clamped( deg, n_cells, mass_line, matrix, row, ind, shift )
929  end do
930  else
931  sll_error('assemble','error in row in assemble3_mass_clamped')
932  end if
933  sll_assert( ind == begin+(2*deg(1)+1)*(2*deg(2)+1)*(2*deg(3)+1))
934 
935  end subroutine assemble_mass3d_clamped
936 
937  !assemble mass matrix in second direction
938  subroutine assemble2_clamped( deg, n_cells, mass_line, matrix, row, ind, shift )
939  sll_int32, intent( in ) :: deg(3)
940  sll_int32, intent( in ) :: n_cells(3)
941  sll_real64, intent( in ) :: mass_line((2*deg(1)+1)*(2*deg(2)+1)*(2*deg(3)+1))
942  type(sll_t_matrix_csr), intent( inout ) :: matrix
943  sll_int32, intent( in ) :: row(3)
944  sll_int32, intent( inout ) :: ind
945  sll_int32, intent( in ) :: shift
946  !local variables
947  sll_int32 :: column, shift1
948  ! For the first deg(2) rows we need to put the first part to the back due to periodic boundaries
949  if(row(2)<=deg(2)) then
950  do column = 2-row(2)+deg(2), deg(2)
951  shift1=(column-1)*(2*deg(1)+1)+shift
952  call assemble1_clamped( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
953  end do
954  do column= deg(2)+1,2*deg(2)+1
955  shift1=(column-1)*(2*deg(1)+1)+shift
956  call assemble1_clamped( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
957  end do
958 
959  do column = 1, deg(2)-row(2)+1
960  shift1=(column-1)*(2*deg(1)+1)+shift
961  call assemble1_clamped( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
962  end do
963  else if(row(2) >= deg(2)+1 .and. row(2)<= n_cells(2)-deg(2)) then
964  do column=1, 2*deg(2)+1
965  shift1=(column-1)*(2*deg(1)+1)+shift
966  call assemble1_clamped( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
967  end do
968  ! For the last deg(2) rows, we need to put the second part to the front due to periodic boundaries
969  else if(row(2) >= n_cells(2)-deg(2)+1 .and. row(2) <= n_cells(2)) then
970 
971  do column = n_cells(2)-row(2)+deg(2)+2, 2*deg(2)+1
972  shift1=(column-1)*(2*deg(1)+1)+shift
973  call assemble1_clamped( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
974  end do
975 
976  do column=1, deg(2)+1
977  shift1=(column-1)*(2*deg(1)+1)+shift
978  call assemble1_clamped(deg, n_cells(1), mass_line, matrix, row(1), ind, shift1)
979  end do
980 
981  do column = 2+deg(2), n_cells(2)-row(2)+deg(2)+1
982  shift1=(column-1)*(2*deg(1)+1)+shift
983  call assemble1_clamped( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
984  end do
985  else
986  sll_error('assemble','error in row in assemble2_mass_clamped')
987  end if
988 
989  end subroutine assemble2_clamped
990 
991  !assemble mass matrix in first direction
992  subroutine assemble1_clamped( deg, n_cells, mass_line, matrix, row, ind, shift )
993  sll_int32, intent( in ) :: deg(3)
994  sll_int32, intent( in ) :: n_cells
995  sll_real64, intent( in ) :: mass_line((2*deg(1)+1)*(2*deg(2)+1)*(2*deg(3)+1))
996  type(sll_t_matrix_csr), intent( inout ) :: matrix
997  sll_int32, intent( in ) :: row
998  sll_int32, intent( inout ) :: ind
999  sll_int32, intent( in ) :: shift
1000 
1001  if(row >= 2*deg(1) .and. row<= n_cells-deg(1)+1) then
1002  matrix%arr_a(1,ind:ind+2*deg(1))=mass_line(1+shift:2*deg(1)+1+shift)
1003 
1004  ind = ind+2*deg(1)+1
1005  else
1006  sll_error('assemble','error in row in assemble1_mass_clamped')
1007  end if
1008 
1009  end subroutine assemble1_clamped
1010 
1012  subroutine assemble_mass3d_clamped_boundary ( deg, n_cells, mass_line, matrix, row, ind )
1013  sll_int32, intent( in ) :: deg(3)
1014  sll_int32, intent( in ) :: n_cells(3)
1015  sll_real64, intent( in ) :: mass_line((7*deg(1)**2-deg(1)-2)/2*(2*deg(2)+1)*(2*deg(3)+1))
1016  type(sll_t_matrix_csr), intent( inout ) :: matrix
1017  sll_int32, intent( in ) :: row(3)
1018  sll_int32, intent( inout ) :: ind
1019  !local variables
1020  sll_int32 :: column, shift
1021 
1022  !assemble massmatrix in third direction
1023  ! For the first deg(3) rows we need to put the first part to the back due to periodic boundaries
1024  if(row(3)<=deg(3)) then
1025  do column = 2-row(3)+deg(3), deg(3)
1026  shift=(column-1)*(7*deg(1)**2-deg(1)-2)/2*(2*deg(2)+1)
1027  call assemble2_clamped_boundary( deg, n_cells, mass_line, matrix, row, ind, shift )
1028  end do
1029  do column= deg(3)+1,2*deg(3)+1
1030  shift=(column-1)*(7*deg(1)**2-deg(1)-2)/2*(2*deg(2)+1)
1031  call assemble2_clamped_boundary( deg, n_cells, mass_line, matrix, row, ind, shift )
1032  end do
1033  do column = 1, deg(3)-row(3)+1
1034  shift=(column-1)*(7*deg(1)**2-deg(1)-2)/2*(2*deg(2)+1)
1035  call assemble2_clamped_boundary( deg, n_cells, mass_line, matrix, row, ind, shift )
1036  end do
1037  else if(row(3) >= deg(3)+1 .and. row(3)<= n_cells(3)-deg(3)) then
1038  do column=1, 2*deg(3)+1
1039  shift=(column-1)*(7*deg(1)**2-deg(1)-2)/2*(2*deg(2)+1)
1040  call assemble2_clamped_boundary( deg, n_cells, mass_line, matrix, row, ind, shift )
1041  end do
1042 
1043  ! For the last deg(3) rows, we need to put the second part to the front due to periodic boundaries
1044  else if(row(3) >= n_cells(3)-deg(3)+1 .and. row(3) <= n_cells(3)) then
1045  do column = n_cells(3)-row(3)+deg(3)+2, 2*deg(3)+1
1046  shift=(column-1)*(7*deg(1)**2-deg(1)-2)/2*(2*deg(2)+1)
1047  call assemble2_clamped_boundary( deg, n_cells, mass_line, matrix, row, ind, shift )
1048  end do
1049  do column=1, deg(3)+1
1050  shift=(column-1)*(7*deg(1)**2-deg(1)-2)/2*(2*deg(2)+1)
1051  call assemble2_clamped_boundary( deg, n_cells, mass_line, matrix, row, ind, shift )
1052  end do
1053  do column = 2+deg(3), n_cells(3)-row(3)+deg(3)+1
1054  shift=(column-1)*(7*deg(1)**2-deg(1)-2)/2*(2*deg(2)+1)
1055  call assemble2_clamped_boundary( deg, n_cells, mass_line, matrix, row, ind, shift )
1056  end do
1057  else
1058  sll_error('assemble','error in row in assemble3_mass_clamped_boundary')
1059  end if
1060 
1061  end subroutine assemble_mass3d_clamped_boundary
1062 
1063  !assemble boundary part in second direction
1064  subroutine assemble2_clamped_boundary( deg, n_cells, mass_line, matrix, row, ind, shift )
1065  sll_int32, intent( in ) :: deg(3)
1066  sll_int32, intent( in ) :: n_cells(3)
1067  sll_real64, intent( in ) :: mass_line((7*deg(1)**2-deg(1)-2)/2*(2*deg(2)+1)*(2*deg(3)+1))
1068  type(sll_t_matrix_csr), intent( inout ) :: matrix
1069  sll_int32, intent( in ) :: row(3)
1070  sll_int32, intent( inout ) :: ind
1071  sll_int32, intent( in ) :: shift
1072  !local variables
1073  sll_int32 :: column, shift1
1074 
1075  ! For the first deg(2) rows we need to put the first part to the back due to periodic boundaries
1076  if(row(2)<=deg(2)) then
1077  do column = 2-row(2)+deg(2), deg(2)
1078  shift1=(column-1)*(7*deg(1)**2-deg(1)-2)/2+shift
1079  call assemble1_clamped_boundary( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1080  end do
1081  do column= deg(2)+1,2*deg(2)+1
1082  shift1=(column-1)*(7*deg(1)**2-deg(1)-2)/2+shift
1083  call assemble1_clamped_boundary( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1084  end do
1085 
1086  do column = 1, deg(2)-row(2)+1
1087  shift1=(column-1)*(7*deg(1)**2-deg(1)-2)/2+shift
1088  call assemble1_clamped_boundary( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1089  end do
1090  else if(row(2) >= deg(2)+1 .and. row(2)<= n_cells(2)-deg(2)) then
1091  do column=1, 2*deg(2)+1
1092  shift1=(column-1)*(7*deg(1)**2-deg(1)-2)/2+shift
1093  call assemble1_clamped_boundary( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1094  end do
1095 
1096  ! For the last deg(2) rows, we need to put the second part to the front due to periodic boundaries
1097  else if(row(2) >= n_cells(2)-deg(2)+1 .and. row(2) <= n_cells(2)) then
1098 
1099  do column = n_cells(2)-row(2)+deg(2)+2, 2*deg(2)+1
1100  shift1=(column-1)*(7*deg(1)**2-deg(1)-2)/2+shift
1101  call assemble1_clamped_boundary( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1102  end do
1103 
1104  do column=1, deg(2)+1
1105  shift1=(column-1)*(7*deg(1)**2-deg(1)-2)/2+shift
1106  call assemble1_clamped_boundary( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1)
1107  end do
1108 
1109  do column = 2+deg(2), n_cells(2)-row(2)+deg(2)+1
1110  shift1=(column-1)*(7*deg(1)**2-deg(1)-2)/2+shift
1111  call assemble1_clamped_boundary( deg, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1112  end do
1113  else
1114  sll_error('assemble','error in row in assemble2_mass_clamped_boundary')
1115  end if
1116 
1117  end subroutine assemble2_clamped_boundary
1118 
1119  !assemble boundary part in first direction
1120  subroutine assemble1_clamped_boundary( deg, n_cells, mass_line, matrix, row, ind, shift )
1121  sll_int32, intent( in ) :: deg(3)
1122  sll_int32, intent( in ) :: n_cells
1123  sll_real64, intent( in ) :: mass_line((7*deg(1)**2-deg(1)-2)/2*(2*deg(2)+1)*(2*deg(3)+1))
1124  type(sll_t_matrix_csr), intent( inout ) :: matrix
1125  sll_int32, intent( in ) :: row
1126  sll_int32, intent( inout ) :: ind
1127  sll_int32, intent( in ) :: shift
1128  !local variables
1129  sll_int32 :: shift1, column
1130 
1131  !use boundary mass for first deg(1) rows
1132  if(row <= deg(1)) then
1133  shift1 = ((row-1)*(2*deg(1)+row))/2 + shift
1134  do column = 1, deg(1)+row
1135  matrix%arr_a(1,ind) = mass_line(column+shift1)
1136  ind = ind+1
1137  end do
1138  else if(row > deg(1) .and. row < 2*deg(1)) then
1139  shift1 = (row-deg(1)-1)*(2*deg(1)+1)+(3*deg(1)**2+deg(1))/2 + shift
1140  do column = 1, 2*deg(1)+1
1141  matrix%arr_a(1,ind) = mass_line(column+shift1)
1142  ind = ind+1
1143  end do
1144  !use boundary mass for last deg(1) rows
1145  else if(row >= n_cells-deg(1)+2 .and. row <= n_cells) then
1146  shift1 = (n_cells-row)*(2*deg(1)+1)+(3*deg(1)**2+deg(1))/2 + shift
1147  do column = 2*deg(1)+1,1,-1
1148  matrix%arr_a(1,ind)=mass_line(column+shift1)
1149  ind = ind+1
1150  end do
1151  else if(row > n_cells .and. row <= n_cells+deg(1)) then
1152  shift1 = ((n_cells+deg(1)-row)*(2*deg(1)+n_cells+deg(1)+1-row))/2 + shift
1153  do column = deg(1)+1+n_cells+deg(1)-row,1,-1
1154  matrix%arr_a(1,ind) = mass_line(column+shift1)
1155  ind = ind+1
1156  end do
1157  else
1158  sll_error('assemble','error in row in assemble1_mass_clamped_boundary')
1159  end if
1160 
1161  end subroutine assemble1_clamped_boundary
1162 
1163 
1164  !---------------------------------------------------------------------------!
1166  subroutine assemble_mixedmass3d ( deg1, deg2, n_cells, mass_line, matrix, row, ind )
1167  sll_int32, intent( in ) :: deg1(3), deg2(3)
1168  sll_int32, intent( in ) :: n_cells(3)
1169  sll_real64, intent( in ) :: mass_line((deg1(1)+deg2(1)+1)*(deg1(2)+ deg2(2)+1)*(deg1(3)+ deg2(3)+1))
1170  type(sll_t_matrix_csr), intent( inout ) :: matrix
1171  sll_int32, intent( in ) :: row(3)
1172  sll_int32, intent( inout ) :: ind
1173  !local variables
1174  sll_int32 :: column, shift
1175 
1176  !assemble massmatrix in third direction
1177  ! For the first deg2(3) rows we need to put the first part to the back due to periodic boundaries
1178  if(row(3) <= deg2(3)) then
1179  do column = 2-row(3)+deg2(3), deg2(3)
1180  shift = (column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1181  call massemble2( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1182  end do
1183  do column = deg2(3)+1,deg1(3)+deg2(3)+1
1184  shift = (column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1185  call massemble2( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1186  end do
1187 
1188  do column = 1, deg2(3)-row(3)+1
1189  shift = (column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1190  call massemble2( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1191  end do
1192 
1193  sll_assert( ind == (row(1)+(row(2)-1)*n_cells(1)+(row(3)-1)*n_cells(1)*n_cells(2))*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)*(deg1(3)+deg2(3)+1)+1 )
1194  else if(row(3) > deg2(3) .and. row(3) <= n_cells(3)-deg1(3)) then
1195  do column = 1, deg1(3)+deg2(3)+1
1196  shift = (column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1197  call massemble2( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1198  end do
1199  sll_assert( ind == (row(1)+(row(2)-1)*n_cells(1)+(row(3)-1)*n_cells(1)*n_cells(2))*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)*(deg1(3)+deg2(3)+1)+1 )
1200  ! For the last deg(3) rows, we need to put the second part to the front due to periodic boundaries
1201  else if(row(3) > n_cells(3)-deg1(3) .and. row(3) <= n_cells(3)) then
1202 
1203  do column = n_cells(3)-row(3)+deg2(3)+2, deg1(3)+deg2(3)+1
1204  shift = (column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1205  call massemble2( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1206  end do
1207 
1208  do column = 1, deg2(3)+1
1209  shift = (column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1210  call massemble2( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1211  end do
1212 
1213  do column = 2+deg2(3), n_cells(3)-row(3)+deg2(3)+1
1214  shift = (column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1215  call massemble2( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1216  end do
1217  sll_assert( ind == (row(1)+(row(2)-1)*n_cells(1)+(row(3)-1)*n_cells(1)*n_cells(2))*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)*(deg1(3)+deg2(3)+1)+1 )
1218  else
1219  sll_error('assemble','error in row in assemble_mass')
1220  end if
1221 
1222  end subroutine assemble_mixedmass3d
1223 
1224  !assemble mass matrix in second direction
1225  subroutine massemble2( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1226  sll_int32, intent( in ) :: deg1(3), deg2(3)
1227  sll_int32, intent( in ) :: n_cells(3)
1228  sll_real64, intent( in ) :: mass_line((deg1(1)+deg2(1)+1)*(deg1(2)+ deg2(2)+1)*(deg1(3)+ deg2(3)+1))
1229  type(sll_t_matrix_csr), intent( inout ) :: matrix
1230  sll_int32, intent( in ) :: row(3)
1231  sll_int32, intent( inout ) :: ind
1232  sll_int32, intent( in ) :: shift
1233  !local variables
1234  sll_int32 :: column, shift1
1235  ! For the first deg2(2) rows we need to put the first part to the back due to periodic boundaries
1236  if(row(2) <= deg2(2)) then
1237  do column = 2-row(2)+deg2(2), deg2(2)
1238  shift1 = (column-1)*(deg1(1)+deg2(1)+1)+shift
1239  call massemble1( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1240  end do
1241  do column = deg2(2)+1,deg1(2)+deg2(2)+1
1242  shift1 = (column-1)*(deg1(1)+deg2(1)+1)+shift
1243  call massemble1( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1244  end do
1245 
1246  do column = 1, deg2(2)-row(2)+1
1247  shift1 = (column-1)*(deg1(1)+deg2(1)+1)+shift
1248  call massemble1( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1249  end do
1250  else if(row(2) > deg2(2) .and. row(2) <= n_cells(2)-deg1(2)) then
1251  do column = 1, deg1(2)+deg2(2)+1
1252  shift1 = (column-1)*(deg1(1)+deg2(1)+1)+shift
1253  call massemble1( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1254  end do
1255 
1256  ! For the last deg1(2) rows, we need to put the second part to the front due to periodic boundaries
1257  else if(row(2) > n_cells(2)-deg1(2) .and. row(2) <= n_cells(2)) then
1258 
1259  do column = n_cells(2)-row(2)+deg2(2)+2, deg1(2)+deg2(2)+1
1260  shift1 = (column-1)*(deg1(1)+deg2(1)+1)+shift
1261  call massemble1( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1262  end do
1263 
1264  do column = 1, deg2(2)+1
1265  shift1 = (column-1)*(deg1(1)+deg2(1)+1)+shift
1266  call massemble1( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1)
1267  end do
1268 
1269  do column = 2+deg2(2), n_cells(2)-row(2)+deg2(2)+1
1270  shift1 = (column-1)*(deg1(1)+deg2(1)+1)+shift
1271  call massemble1( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1272  end do
1273  else
1274  sll_error('assemble','error in row in assemble_mass')
1275  end if
1276 
1277  end subroutine massemble2
1278 
1279  !assemble mass matrix in first direction
1280  subroutine massemble1( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1281  sll_int32, intent( in ) :: deg1(3), deg2(3)
1282  sll_int32, intent( in ) :: n_cells
1283  sll_real64, intent( in ) :: mass_line((deg1(1)+deg2(1)+1)*(deg1(2)+ deg2(2)+1)*(deg1(3)+ deg2(3)+1))
1284  type(sll_t_matrix_csr), intent( inout ) :: matrix
1285  sll_int32, intent( in ) :: row
1286  sll_int32, intent( inout ) :: ind
1287  sll_int32, intent( in ) :: shift
1288  !local variables
1289  sll_int32 :: column
1290 
1291  ! For the first deg2(1) rows we need to put the first part to the back due to periodic boundaries
1292  if(row <= deg2(1)) then
1293  do column = 2-row+deg2(1), deg2(1)
1294  matrix%arr_a(1,ind) = mass_line(column+shift)
1295  ind = ind+1
1296  end do
1297  matrix%arr_a(1,ind:ind+deg1(1)) = mass_line(deg2(1)+1+shift:deg1(1)+deg2(1)+1+shift)
1298  ind = ind+deg1(1)+1
1299 
1300  do column = 1, deg2(1)-row+1
1301  matrix%arr_a(1,ind) = mass_line(column+shift)
1302  ind = ind+1
1303  end do
1304  else if(row > deg2(1) .and. row <= n_cells-deg1(1)) then
1305  matrix%arr_a(1,ind:ind+deg1(1)+deg2(1))=mass_line(1+shift:deg1(1)+deg2(1)+1+shift)
1306 
1307  ind = ind+deg1(1)+deg2(1)+1
1308  ! For the last deg1(1) rows, we need to put the second part to the front due to periodic boundaries
1309  else if(row > n_cells-deg1(1) .and. row <= n_cells) then
1310 
1311  do column = n_cells-row+deg2(1)+2, deg1(1)+deg2(1)+1
1312  matrix%arr_a(1,ind) = mass_line(column+shift)
1313  ind = ind+1
1314  end do
1315 
1316  matrix%arr_a(1,ind:ind+deg2(1)) = mass_line(1+shift:deg2(1)+1+shift)
1317  ind = ind+deg2(1)+1
1318 
1319  do column = 2+deg2(1), n_cells-row+deg2(1)+1
1320  matrix%arr_a(1,ind) = mass_line(column+shift)
1321  ind = ind+1
1322  end do
1323  else
1324  sll_error('assemble','error in row in assemble_mass')
1325  end if
1326 
1327  end subroutine massemble1
1328 
1329  !---------------------------------------------------------------------------!
1331  subroutine assemble_mixedmass3d_clamped ( deg1, deg2, n_cells, mass_line, matrix, row, ind )
1332  sll_int32, intent( in ) :: deg1(3), deg2(3)
1333  sll_int32, intent( in ) :: n_cells(3)
1334  sll_real64, intent( in ) :: mass_line((deg1(1)+deg2(1)+1)*(deg1(2)+ deg2(2)+1)*(deg1(3)+ deg2(3)+1))
1335  type(sll_t_matrix_csr), intent( inout ) :: matrix
1336  sll_int32, intent( in ) :: row(3)
1337  sll_int32, intent( inout ) :: ind
1338  !local variables
1339  sll_int32 :: column, shift, begin
1340 
1341  begin= ind
1342 
1343  !assemble massmatrix in third direction
1344  ! For the first deg2(3) rows we need to put the first part to the back due to periodic boundaries
1345  if(row(3)<=deg2(3)) then
1346  do column = 2-row(3)+deg2(3), deg2(3)
1347  shift=(column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1348  call massemble2_clamped( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1349  end do
1350  do column= deg2(3)+1,deg1(3)+deg2(3)+1
1351  shift=(column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1352  call massemble2_clamped( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1353  end do
1354 
1355  do column = 1, deg2(3)-row(3)+1
1356  shift=(column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1357  call massemble2_clamped( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1358  end do
1359  else if(row(3) > deg2(3) .and. row(3)<= n_cells(3)-deg1(3)) then
1360  do column=1, deg1(3)+deg2(3)+1
1361  shift=(column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1362  call massemble2_clamped( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1363  end do
1364  ! For the last deg(3) rows, we need to put the second part to the front due to periodic boundaries
1365  else if(row(3) > n_cells(3)-deg1(3) .and. row(3) <= n_cells(3)) then
1366 
1367  do column = n_cells(3)-row(3)+deg2(3)+2, deg1(3)+deg2(3)+1
1368  shift=(column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1369  call massemble2_clamped( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1370  end do
1371 
1372  do column=1, deg2(3)+1
1373  shift=(column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1374  call massemble2_clamped( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1375  end do
1376 
1377  do column = 2+deg2(3), n_cells(3)-row(3)+deg2(3)+1
1378  shift=(column-1)*(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)
1379  call massemble2_clamped( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1380  end do
1381  else
1382  sll_error('assemble','error in row in assemble_mass')
1383  end if
1384  sll_assert( ind == begin+(deg1(1)+deg2(1)+1)*(deg1(2)+deg2(2)+1)*(deg1(3)+deg2(3)+1) )
1385 
1386  end subroutine assemble_mixedmass3d_clamped
1387 
1388  !assemble mass matrix in second direction
1389  subroutine massemble2_clamped( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1390  sll_int32, intent( in ) :: deg1(3), deg2(3)
1391  sll_int32, intent( in ) :: n_cells(3)
1392  sll_real64, intent( in ) :: mass_line((deg1(1)+deg2(1)+1)*(deg1(2)+ deg2(2)+1)*(deg1(3)+ deg2(3)+1))
1393  type(sll_t_matrix_csr), intent( inout ) :: matrix
1394  sll_int32, intent( in ) :: row(3)
1395  sll_int32, intent( inout ) :: ind
1396  sll_int32, intent( in ) :: shift
1397  !local variables
1398  sll_int32 :: column, shift1
1399 
1400  ! For the first deg2(2) rows we need to put the first part to the back due to periodic boundaries
1401  if(row(2)<=deg2(2)) then
1402  do column = 2-row(2)+deg2(2), deg2(2)
1403  shift1=(column-1)*(deg1(1)+deg2(1)+1)+shift
1404  call massemble1_clamped( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1405  end do
1406  do column= deg2(2)+1,deg1(2)+deg2(2)+1
1407  shift1=(column-1)*(deg1(1)+deg2(1)+1)+shift
1408  call massemble1_clamped( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1409  end do
1410 
1411  do column = 1, deg2(2)-row(2)+1
1412  shift1=(column-1)*(deg1(1)+deg2(1)+1)+shift
1413  call massemble1_clamped( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1414  end do
1415  else if(row(2) > deg2(2) .and. row(2)<= n_cells(2)-deg1(2)) then
1416  do column=1, deg1(2)+deg2(2)+1
1417  shift1=(column-1)*(deg1(1)+deg2(1)+1)+shift
1418  call massemble1_clamped( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1419  end do
1420  ! For the last deg1(2) rows, we need to put the second part to the front due to periodic boundaries
1421  else if(row(2) > n_cells(2)-deg1(2) .and. row(2) <= n_cells(2)) then
1422 
1423  do column = n_cells(2)-row(2)+deg2(2)+2, deg1(2)+deg2(2)+1
1424  shift1=(column-1)*(deg1(1)+deg2(1)+1)+shift
1425  call massemble1_clamped( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1426  end do
1427 
1428  do column=1, deg2(2)+1
1429  shift1=(column-1)*(deg1(1)+deg2(1)+1)+shift
1430  call massemble1_clamped( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1431  end do
1432 
1433  do column = 2+deg2(2), n_cells(2)-row(2)+deg2(2)+1
1434  shift1=(column-1)*(deg1(1)+deg2(1)+1)+shift
1435  call massemble1_clamped( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1436  end do
1437  else
1438  sll_error('assemble','error in row in assemble_mass')
1439  end if
1440 
1441  end subroutine massemble2_clamped
1442 
1443  !assemble mass matrix in first direction
1444  subroutine massemble1_clamped( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1445  sll_int32, intent( in ) :: deg1(3), deg2(3)
1446  sll_int32, intent( in ) :: n_cells
1447  sll_real64, intent( in ) :: mass_line((deg1(1)+deg2(1)+1)*(deg1(2)+ deg2(2)+1)*(deg1(3)+ deg2(3)+1))
1448  type(sll_t_matrix_csr), intent( inout ) :: matrix
1449  sll_int32, intent( in ) :: row
1450  sll_int32, intent( inout ) :: ind
1451  sll_int32, intent( in ) :: shift
1452 
1453  if(row > deg1(1)+deg2(1) .and. row <= n_cells-deg2(1)) then
1454  matrix%arr_a(1,ind:ind+deg1(1)+deg2(1))=mass_line(1+shift:deg1(1)+deg2(1)+1+shift)
1455  ind = ind+deg1(1)+deg2(1)+1
1456  else
1457  sll_error('assemble','error in row in assemble1_mixedmass_clamped')
1458  end if
1459 
1460  end subroutine massemble1_clamped
1461 
1462  !---------------------------------------------------------------------------!
1464  subroutine assemble_mixedmass3d_clamped_boundary ( deg1, deg2, n_cells, mass_line, matrix, row, ind )
1465  sll_int32, intent( in ) :: deg1(3), deg2(3)
1466  sll_int32, intent( in ) :: n_cells(3)
1467  sll_real64, intent( in ) :: mass_line(((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)*(deg1(2)+deg2(2)+1)*(deg1(3)+deg2(3)+1))
1468  type(sll_t_matrix_csr), intent( inout ) :: matrix
1469  sll_int32, intent( in ) :: row(3)
1470  sll_int32, intent( inout ) :: ind
1471  !local variables
1472  sll_int32 :: column, shift
1473 
1474  !assemble massmatrix in third direction
1475  ! For the first deg2(3) rows we need to put the first part to the back due to periodic boundaries
1476  if(row(3)<=deg2(3)) then
1477  do column = 2-row(3)+deg2(3), deg2(3)
1478  shift=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)*(deg1(2)+deg2(2)+1)
1479  call massemble2_clamped_boundary( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1480  end do
1481  do column= deg2(3)+1,deg1(3)+deg2(3)+1
1482  shift=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)*(deg1(2)+deg2(2)+1)
1483  call massemble2_clamped_boundary( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1484  end do
1485 
1486  do column = 1, deg2(3)-row(3)+1
1487  shift=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)*(deg1(2)+deg2(2)+1)
1488  call massemble2_clamped_boundary( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1489  end do
1490  else if(row(3) > deg2(3) .and. row(3)<= n_cells(3)-deg1(3)) then
1491  do column=1, deg1(3)+deg2(3)+1
1492  shift=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)*(deg1(2)+deg2(2)+1)
1493  call massemble2_clamped_boundary( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1494  end do
1495  ! For the last deg(3) rows, we need to put the second part to the front due to periodic boundaries
1496  else if(row(3) > n_cells(3)-deg1(3) .and. row(3) <= n_cells(3)) then
1497 
1498  do column = n_cells(3)-row(3)+deg2(3)+2, deg1(3)+deg2(3)+1
1499  shift=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)*(deg1(2)+deg2(2)+1)
1500  call massemble2_clamped_boundary( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1501  end do
1502 
1503  do column=1, deg2(3)+1
1504  shift=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)*(deg1(2)+deg2(2)+1)
1505  call massemble2_clamped_boundary( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1506  end do
1507 
1508  do column = 2+deg2(3), n_cells(3)-row(3)+deg2(3)+1
1509  shift=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)*(deg1(2)+deg2(2)+1)
1510  call massemble2_clamped_boundary( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1511  end do
1512  else
1513  sll_error('assemble','error in row in assemble_mass')
1514  end if
1515 
1516  end subroutine assemble_mixedmass3d_clamped_boundary
1517 
1518  !assemble boundary part in second direction
1519  subroutine massemble2_clamped_boundary( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1520  sll_int32, intent( in ) :: deg1(3), deg2(3)
1521  sll_int32, intent( in ) :: n_cells(3)
1522  sll_real64, intent( in ) :: mass_line(((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)*(deg1(2)+ deg2(2)+1)*(deg1(3)+ deg2(3)+1))
1523  type(sll_t_matrix_csr), intent( inout ) :: matrix
1524  sll_int32, intent( in ) :: row(3)
1525  sll_int32, intent( inout ) :: ind
1526  sll_int32, intent( in ) :: shift
1527  !local variables
1528  sll_int32 :: column, shift1
1529 
1530  ! For the first deg2(2) rows we need to put the first part to the back due to periodic boundaries
1531  if(row(2)<=deg2(2)) then
1532  do column = 2-row(2)+deg2(2), deg2(2)
1533  shift1=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)+shift
1534  call massemble1_clamped_boundary( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1535  end do
1536  do column= deg2(2)+1,deg1(2)+deg2(2)+1
1537  shift1=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)+shift
1538  call massemble1_clamped_boundary( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1539  end do
1540 
1541  do column = 1, deg2(2)-row(2)+1
1542  shift1=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)+shift
1543  call massemble1_clamped_boundary( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1544  end do
1545  else if(row(2) > deg2(2) .and. row(2)<= n_cells(2)-deg1(2)) then
1546  do column=1, deg1(2)+deg2(2)+1
1547  shift1=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)+shift
1548  call massemble1_clamped_boundary( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1549  end do
1550  ! For the last deg1(2) rows, we need to put the second part to the front due to periodic boundaries
1551  else if(row(2) > n_cells(2)-deg1(2) .and. row(2) <= n_cells(2)) then
1552 
1553  do column = n_cells(2)-row(2)+deg2(2)+2, deg1(2)+deg2(2)+1
1554  shift1=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)+shift
1555  call massemble1_clamped_boundary( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1556  end do
1557 
1558  do column=1, deg2(2)+1
1559  shift1=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)+shift
1560  call massemble1_clamped_boundary( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1)
1561  end do
1562 
1563  do column = 2+deg2(2), n_cells(2)-row(2)+deg2(2)+1
1564  shift1=(column-1)*((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)+shift
1565  call massemble1_clamped_boundary( deg1, deg2, n_cells(1), mass_line, matrix, row(1), ind, shift1 )
1566  end do
1567  else
1568  sll_error('assemble','error in row in assemble_mass')
1569  end if
1570 
1571  end subroutine massemble2_clamped_boundary
1572 
1573  !assemble boundary part in first direction
1574  subroutine massemble1_clamped_boundary( deg1, deg2, n_cells, mass_line, matrix, row, ind, shift )
1575  sll_int32, intent( in ) :: deg1(3), deg2(3)
1576  sll_int32, intent( in ) :: n_cells
1577  sll_real64, intent( in ) :: mass_line(((deg1(1)+deg2(1))**2+(2*deg2(1)+deg1(1)-deg1(1)**2)/2)*(deg1(2)+ deg2(2)+1)*(deg1(3)+ deg2(3)+1))
1578  type(sll_t_matrix_csr), intent( inout ) :: matrix
1579  sll_int32, intent( in ) :: row
1580  sll_int32, intent( inout ) :: ind
1581  sll_int32, intent( in ) :: shift
1582  !local variables
1583  sll_int32 :: column, shift1
1584 
1585  !use boundary mass for first deg1(1) rows
1586  if(row<=deg1(1)) then
1587  shift1 = ((row-1)*(2*deg2(1)+row))/2 +shift
1588  do column = 1, deg2(1)+row
1589  matrix%arr_a(1,ind) = mass_line(column+shift1)
1590  ind = ind+1
1591  end do
1592  else if(row > deg1(1) .and. row <= deg1(1)+deg2(1)) then
1593  shift1 = (row-deg1(1)-1)*(deg1(1)+deg2(1)+1)+(deg1(1)*(2*deg2(1)+deg1(1)+1))/2 + shift
1594  do column = 1, deg1(1)+deg2(1)+1
1595  matrix%arr_a(1,ind) = mass_line(column+shift1)
1596  ind = ind+1
1597  end do
1598 
1599  !use boundary mass for last deg1(1) rows
1600  else if(row >= n_cells-deg2(1)+1 .and. row <= n_cells) then
1601  shift1 = (n_cells-row)*(deg1(1)+deg2(1)+1)+(deg1(1)*(2*deg2(1)+deg1(1)+1))/2 + shift
1602  do column = deg1(1)+deg2(1)+1,1,-1
1603  matrix%arr_a(1,ind)=mass_line(column+shift1)
1604  ind = ind+1
1605  end do
1606  else if(row > n_cells .and. row <= n_cells+deg1(1)) then
1607  shift1 = ((n_cells+deg1(1)-row)*(2*deg2(1)+n_cells+deg1(1)+1-row))/2 + shift
1608  do column = deg2(1)+1+n_cells+deg1(1)-row,1,-1
1609  matrix%arr_a(1,ind) = mass_line(column+shift1)
1610  ind = ind+1
1611  end do
1612  else
1613  sll_error('assemble','error in row in assemble_mass')
1614  end if
1615 
1616  end subroutine massemble1_clamped_boundary
1617 
1618 
1619 end module sll_m_spline_fem_utilities_3d_helper
sll_m_matrix_csr
module for Compressed Sparse Row Matrix (CSR)
Definition: sll_m_matrix_csr.F90:12
sll_m_spline_fem_utilities_3d_helper
Helper for spline finite elements utilites.
Definition: sll_m_spline_fem_utilities_3d_helper.F90:8
sll_m_spline_fem_utilities_3d_helper::mloop2
subroutine mloop2(i, j, shift, deg1, deg2, n_cells, ind, spmat)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:461
sll_m_spline_fem_utilities_3d_helper::massemble2_clamped_boundary
subroutine massemble2_clamped_boundary(deg1, deg2, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1520
sll_m_spline_fem_utilities_3d_helper::assemble1_clamped
subroutine assemble1_clamped(deg, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:993
sll_m_spline_fem_utilities_3d_helper::massemble1_clamped
subroutine massemble1_clamped(deg1, deg2, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1445
sll_m_spline_fem_utilities_3d_helper::loop2
subroutine loop2(i, j, shift, deg, n_cells, ind, spmat)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:120
sll_m_spline_fem_utilities_3d_helper::assemble_mixedmass3d
subroutine, public assemble_mixedmass3d(deg1, deg2, n_cells, mass_line, matrix, row, ind)
Assemble the given row of the mixed mass matrix.
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1167
sll_m_spline_fem_utilities_3d_helper::assemble_mixedmass3d_clamped
subroutine, public assemble_mixedmass3d_clamped(deg1, deg2, n_cells, mass_line, matrix, row, ind)
Assemble the given row of the clamped mixed mass matrix.
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1332
sll_m_spline_fem_utilities_3d_helper::mloop1
subroutine mloop1(i, shift, deg1, deg2, n_cells, ind, spmat)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:505
sll_m_spline_fem_utilities_3d_helper::massemble2_clamped
subroutine massemble2_clamped(deg1, deg2, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1390
sll_m_spline_fem_utilities_3d_helper::assemble2_clamped
subroutine assemble2_clamped(deg, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:939
sll_m_spline_fem_utilities_3d_helper::loop1
subroutine loop1(i, shift, deg, n_cells, ind, spmat)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:164
sll_m_spline_fem_utilities_3d_helper::assemble2_clamped_boundary
subroutine assemble2_clamped_boundary(deg, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1065
sll_m_spline_fem_utilities_3d_helper::loop1_clamped
subroutine loop1_clamped(i, shift, deg, n_cells, ind, spmat)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:346
sll_m_spline_fem_utilities_3d_helper::assemble_mass3d
subroutine, public assemble_mass3d(deg, n_cells, mass_line, matrix, row, ind)
Assemble the given row of the 3d mass matrix.
Definition: sll_m_spline_fem_utilities_3d_helper.F90:715
sll_m_spline_fem_utilities_3d_helper::sll_s_spline_fem_sparsity_mixedmass3d_clamped
subroutine, public sll_s_spline_fem_sparsity_mixedmass3d_clamped(deg1, deg2, n_cells, spmat)
Helper function to create sparsity pattern of the 3d clamped mixed mass matrix.
Definition: sll_m_spline_fem_utilities_3d_helper.F90:548
sll_m_spline_fem_utilities_3d_helper::massemble1_clamped_boundary
subroutine massemble1_clamped_boundary(deg1, deg2, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1575
sll_m_spline_fem_utilities_3d_helper::assemble_mass3d_clamped_boundary
subroutine, public assemble_mass3d_clamped_boundary(deg, n_cells, mass_line, matrix, row, ind)
Assemble the boundary part of the clamped mass matrix.
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1013
sll_m_spline_fem_utilities_3d_helper::sll_s_spline_fem_sparsity_mass3d
subroutine, public sll_s_spline_fem_sparsity_mass3d(deg, n_cells, spmat)
Helper function to create sparsity pattern of the 3d mass matrix.
Definition: sll_m_spline_fem_utilities_3d_helper.F90:40
sll_m_spline_fem_utilities_3d_helper::assemble1
subroutine assemble1(deg, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:828
sll_m_spline_fem_utilities_3d_helper::loop2_clamped
subroutine loop2_clamped(i, j, shift, deg, n_cells, ind, spmat)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:299
sll_m_spline_fem_utilities_3d_helper::massemble2
subroutine massemble2(deg1, deg2, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1226
sll_m_spline_fem_utilities_3d_helper::assemble_mass3d_clamped
subroutine, public assemble_mass3d_clamped(deg, n_cells, mass_line, matrix, row, ind)
Assemble the given row of the clamped mass matrix.
Definition: sll_m_spline_fem_utilities_3d_helper.F90:880
sll_m_spline_fem_utilities_3d_helper::assemble1_clamped_boundary
subroutine assemble1_clamped_boundary(deg, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1121
sll_m_spline_fem_utilities_3d_helper::assemble_mixedmass3d_clamped_boundary
subroutine, public assemble_mixedmass3d_clamped_boundary(deg1, deg2, n_cells, mass_line, matrix, row, ind)
Assemble the boundary part of the clamped mixed mass matrix.
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1465
sll_m_spline_fem_utilities_3d_helper::mloop1_clamped
subroutine mloop1_clamped(i, shift, deg1, deg2, n_cells, ind, spmat)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:681
sll_m_spline_fem_utilities_3d_helper::sll_s_spline_fem_sparsity_mixedmass3d
subroutine, public sll_s_spline_fem_sparsity_mixedmass3d(deg1, deg2, n_cells, spmat)
Helper function to create sparsity pattern of the 3d mass matrix.
Definition: sll_m_spline_fem_utilities_3d_helper.F90:383
sll_m_spline_fem_utilities_3d_helper::sll_s_spline_fem_sparsity_mass3d_clamped
subroutine, public sll_s_spline_fem_sparsity_mass3d_clamped(deg, n_cells, spmat)
Helper function to create sparsity pattern of the 3d clamped mass matrix.
Definition: sll_m_spline_fem_utilities_3d_helper.F90:207
sll_m_spline_fem_utilities_3d_helper::massemble1
subroutine massemble1(deg1, deg2, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:1281
sll_m_spline_fem_utilities_3d_helper::mloop2_clamped
subroutine mloop2_clamped(i, j, shift, deg1, deg2, n_cells, ind, spmat)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:637
sll_m_spline_fem_utilities_3d_helper::assemble2
subroutine assemble2(deg, n_cells, mass_line, matrix, row, ind, shift)
Definition: sll_m_spline_fem_utilities_3d_helper.F90:773
sll_m_matrix_csr::sll_t_matrix_csr
class for a csr matrix
Definition: sll_m_matrix_csr.F90:33
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1