GCC Code Coverage Report


Directory: src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 97.8% 452 / 0 / 462
Functions: 100.0% 22 / 0 / 22
Branches: 62.5% 565 / 0 / 904

mform/m_mform_constraint_polar.f90
Line Branch Exec Source
1 !> @brief Module containing implementations of various local m-form constraints, such as:
2 !> - MFormPolarRegularity, which imposes a polar regularity constraint on an m-form
3 module m_mform_constraint_polar
4 use m_bspline, only: BSplineFun
5 use m_common, only: wp
6 use m_mform_basis, only: MFormSpace
7 use m_sparsemat, only: SparseMat
8 use m_tensorprod_indices, only: TensorProdIndexList, size, size_2d
9 use m_tensorprod_shared, only: SharedMemoryWindow
10 use m_mform_constraint_abstract, only: MFormConstraintLocal
11 #include "petsc.fi"
12
13 implicit none
14
15 private
16 public :: MFormPolarRegularity, MFormPolarBoundedness
17
18 !> The threshold under which extraction/projection coefficients are considered zero
19 real(wp), parameter :: COEFF_THR = 1e-14_wp
20
21 !> An object that represents a polar regularity constraint on an m-form space
22 type, extends(MFormConstraintLocal) :: MFormPolarRegularity
23 !> The regularity of the polar constraint
24 integer :: regularity
25 !> Enriched space
26 logical :: enriched
27 !> Orthogonal extraction (Euclidean inner product)
28 logical :: orthogonal
29 !> The radial B-splines for the polar regularity constraint
30 type(BSplineFun), allocatable :: radial_bsplines(:)
31 !> The angular B-splines for the polar regularity constraint
32 type(BSplineFun), allocatable :: angular_bsplines(:)
33
34 integer, private :: row0_2d(4)
35
36 !> The zero-form space
37 type(MFormSpace) :: space0
38
39 !> Communicator for the shared memory window
40 integer :: comm_window = -1, my_rank_window = -1
41
42 !> The dense 2D extractor operators (shared memory)
43 real(wp), pointer :: E1(:, :, :, :), E2(:, :, :, :) ! i, j, block, polar_index (0, 0, 1, 1)-based indexing
44
45 !> The shared memory windows for the extractor matrices
46 type(SharedMemoryWindow) :: shmem_E1, shmem_E2
47 contains
48 procedure :: init => init_mform_polar
49 procedure :: apply => apply_mform_polar
50 procedure :: destroy_derived => destroy_mform_polarregularity
51 procedure, private :: init_shmem_window
52 procedure, private :: compute_1d_bsplines
53 procedure, private :: insert_polar_extractor_0forms
54 procedure, private :: insert_polar_extractor_1forms
55 procedure, private :: insert_polar_extractor_1forms_enrich
56 procedure, private :: insert_polar_extractor_2forms
57 procedure, private :: insert_polar_extractor_2forms_enrich
58 procedure, private :: insert_polar_extractor_3forms
59 procedure, private :: insert_polar_extractor_3forms_enrich
60 procedure, private :: finalize_polar_extractor_0forms
61 procedure, private :: finalize_polar_extractor_1forms
62 procedure, private :: finalize_polar_extractor_2forms
63 procedure, private :: finalize_polar_extractor_3forms
64 end type
65
66 interface MFormPolarRegularity
67 module procedure construct_mform_polar
68 end interface
69
70 !> An alias for MFormPolarRegularity(0), i.e., a boundedness constraint on an m-form resulting in conforming FEM spaces
71 interface MFormPolarBoundedness
72 module procedure construct_mform_bounded
73 end interface
74 contains
75
76 !> @brief Construct a boundedness constraint on an m-form
77 !>
78 !> @return The constructed MFormPolarRegularity object
79 40 pure type(MFormPolarRegularity) function construct_mform_bounded() result(this)
80 implicit none
81
82 this = MFormPolarRegularity(0)
83 40 end function construct_mform_bounded
84
85 !> @brief Construct a polar regularity constraint on an m-form
86 !>
87 !> @param[in] _(optional)_ regularity The regularity of the polar constraint (default: degree of the B-spline basis in the
88 !> x-direction)
89 !> @param[in] _(optional)_ enriched Whether to use an enriched polar regularity constraint (default: .true.)
90 !> @param[in] _(optional)_ orthogonal Whether to use orthogonal extraction operators (default: .true.)
91 !>
92 !> @return The constructed MFormPolarRegularity object
93 946 pure type(MFormPolarRegularity) function construct_mform_polar(regularity, enriched, orthogonal) result(this)
94 implicit none
95
96 integer, intent(in), optional :: regularity
97 logical, intent(in), optional :: enriched
98 logical, intent(in), optional :: orthogonal
99
100
2/2
✓ Branch 2 → 3 taken 834 times.
✓ Branch 2 → 4 taken 152 times.
986 if (present(regularity)) then
101 794 this%regularity = regularity
102 else
103 this%regularity = -1 ! will be set properly in init method based on degree in x-direction
104 end if
105
106
2/2
✓ Branch 4 → 5 taken 818 times.
✓ Branch 4 → 6 taken 128 times.
946 if (present(enriched)) then
107 818 this%enriched = enriched
108 else
109 this%enriched = .true.
110 end if
111
112
2/2
✓ Branch 6 → 7 taken 280 times.
✓ Branch 6 → 8 taken 666 times.
946 if (present(orthogonal)) then
113 280 this%orthogonal = orthogonal
114 else
115 this%orthogonal = .true.
116 end if
117 946 end function construct_mform_polar
118
119 !> @brief Initialize a polar regularity constraint on an m-form space
120 !>
121 !> @param[inout] this The MFormPolarRegularity object to initialize
122 !> @param[in] space The MFormSpace to which the polar regularity constraint applies
123 !> @param[in] _(optional)_ coord_transform The coordinate transformation of the m-form space
124 1321 subroutine init_mform_polar(this, space, coord_transform)
125 use m_mform_basis, only: MFormSpace, get_zero_form_space
126 use m_qr_factorisation, only: orthogonalize
127 use m_tensorprod, only: size, TensorProdIndices
128 use m_coord_transform_abstract, only: CoordTransformAbstract
129 implicit none
130
131 class(MFormPolarRegularity), intent(inout) :: this
132 type(MFormSpace), intent(in) :: space
133 class(CoordTransformAbstract), optional, intent(in) :: coord_transform
134
135 integer :: s1, s2, nr_blocks
136
137 logical :: coord_has_polar_xy_singularity_
138
139 type(TensorProdIndices) :: inds
140
141 1321 call this%destroy()
142
143
5/8
✓ Branch 3 → 4 taken 1321 times.
✗ Branch 3 → 7 not taken.
✓ Branch 4 → 5 taken 1321 times.
✗ Branch 4 → 6 not taken.
✓ Branch 7 → 8 taken 1321 times.
✗ Branch 7 → 10 not taken.
✓ Branch 8 → 9 taken 132 times.
✓ Branch 8 → 10 taken 1189 times.
1321 this%space = space
144
2/2
✓ Branch 11 → 12 taken 132 times.
✓ Branch 11 → 13 taken 1189 times.
1321 this%space0 = get_zero_form_space(space)
145
146 ! The constraint is shared amongst MPI ranks, however, most ranks contribute zeros to the extractor
147 1321 this%is_mpi_shared(1) = this%space%tp_spaces(1)%domain%nr_subintervals(1) > 1
148 1321 this%is_mpi_shared(2) = this%space%tp_spaces(1)%domain%nr_subintervals(2) > 1
149 1321 this%is_mpi_shared(3) = .false.
150
151 coord_has_polar_xy_singularity_ = .false.
152
3/4
✓ Branch 13 → 14 taken 1277 times.
✓ Branch 13 → 16 taken 44 times.
✓ Branch 14 → 15 taken 1277 times.
✗ Branch 14 → 16 not taken.
1321 if (present(coord_transform)) coord_has_polar_xy_singularity_ = coord_transform%has_polar_xy_singularity
153
154 1321 nr_blocks = size(space%tp_spaces)
155
9/16
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 1321 times.
✓ Branch 18 → 17 taken 1321 times.
✗ Branch 18 → 19 not taken.
✓ Branch 19 → 20 taken 1321 times.
✗ Branch 19 → 21 not taken.
✓ Branch 21 → 22 taken 1321 times.
✗ Branch 21 → 23 not taken.
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 1321 times.
✗ Branch 25 → 26 not taken.
✓ Branch 25 → 27 taken 1321 times.
✗ Branch 27 → 28 not taken.
✓ Branch 27 → 29 taken 1321 times.
✓ Branch 30 → 31 taken 2561 times.
✓ Branch 30 → 32 taken 1321 times.
6524 allocate (this%index_list(nr_blocks))
156
9/16
✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 1321 times.
✓ Branch 34 → 33 taken 1321 times.
✗ Branch 34 → 35 not taken.
✓ Branch 35 → 36 taken 1321 times.
✗ Branch 35 → 37 not taken.
✓ Branch 37 → 38 taken 1321 times.
✗ Branch 37 → 39 not taken.
✗ Branch 39 → 40 not taken.
✓ Branch 39 → 41 taken 1321 times.
✗ Branch 41 → 42 not taken.
✓ Branch 41 → 43 taken 1321 times.
✗ Branch 43 → 44 not taken.
✓ Branch 43 → 45 taken 1321 times.
✓ Branch 46 → 47 taken 2561 times.
✓ Branch 46 → 48 taken 1321 times.
6524 allocate (this%extractor(nr_blocks))
157
158
2/2
✓ Branch 48 → 49 taken 44 times.
✓ Branch 48 → 58 taken 1277 times.
1321 if (.not. coord_has_polar_xy_singularity_) then
159 44 call inds%empty()
160
161
2/2
✓ Branch 51 → 52 taken 44 times.
✓ Branch 51 → 56 taken 44 times.
88 do s1 = 1, nr_blocks
162 44 call this%index_list(s1)%init(inds, space%tp_spaces(s1)%rank_resp_bspline, space%tp_spaces(s1)%rank_l0)
163 44 call this%extractor(s1)%init(nr_rows=0, nr_cols=0, nr_nonzero=0)
164 88 call this%extractor(s1)%finalize()
165 end do
166
167 44 return
168 end if
169
170
1/4
✗ Branch 58 → 59 not taken.
✓ Branch 58 → 65 taken 1277 times.
✗ Branch 59 → 60 not taken.
✗ Branch 59 → 65 not taken.
1277 if (space%tp_spaces(1)%domain%is_distributed_memory(1) .and. space%tp_spaces(1)%domain%my_subinterval_ijk(1) > 0) then
171 do s1 = 1, size(space%tp_spaces)
172 if (space%tp_spaces(s1)%rank_resp_bspline%i1 < this%regularity) then
173 error stop "MFormPolarRegularity::init_mform_polar: the first subinterval in the x-direction of the tensor product"// &
174 & " space must have at least 'regularity + 1' B-splines whenever this direction is distributed in memory."
175 end if
176 end do
177 end if
178
1/2
✗ Branch 65 → 66 not taken.
✓ Branch 65 → 67 taken 1277 times.
1277 if (space%tp_spaces(1)%domain%is_distributed_memory(2)) then
179 error stop "MFormPolarRegularity::init_mform_polar: The polar regularity constraint is not compatible with a distributed "// &
180 & "memory domain decomposition in the y-direction"
181 end if
182
183
2/4
✓ Branch 67 → 68 taken 1277 times.
✗ Branch 67 → 69 not taken.
✗ Branch 68 → 69 not taken.
✓ Branch 68 → 70 taken 1277 times.
1277 if (space%tp_spaces(1)%spaces(1)%is_periodic .or. .not. space%tp_spaces(1)%spaces(2)%is_periodic) then
184 error stop "MFormPolarRegularity::init_mform_polar: The polar regularity constraint requires the B-spline basis in the"// &
185 & " x-direction to be non-periodic and the B-spline basis in the y-direction to be periodic"
186 end if
187
188 1277 this%is_homogeneous = .true.
189 1277 this%is_symmetric = .false.
190
191
2/2
✓ Branch 70 → 71 taken 160 times.
✓ Branch 70 → 72 taken 1117 times.
1277 if (this%regularity <= -1) then
192 160 this%regularity = this%space0%tp_spaces(1)%spaces(1)%degree
193 else
194
1/2
✗ Branch 72 → 73 not taken.
✓ Branch 72 → 74 taken 1117 times.
1117 if (this%regularity < 0 .or. this%regularity > this%space0%tp_spaces(1)%spaces(1)%degree) then
195 error stop "MFormPolarRegularity::init_mform_polar: 'regularity' must be in [0, degree_x] where 'degree_x' is the "// &
196 & "degree of the B-spline basis in the x-direction of the zero-form space"
197 end if
198 end if
199
200 ! The forward part of commutativty (strong part of Hodge decomposition)
201 1810 select case (space%m)
202 case (0)
203 533 call this%insert_polar_extractor_0forms()
204 case (1)
205 371 call this%insert_polar_extractor_1forms()
206 case (2)
207 249 call this%insert_polar_extractor_2forms()
208 case (3)
209
4/5
✓ Branch 74 → 75 taken 533 times.
✓ Branch 74 → 76 taken 371 times.
✓ Branch 74 → 77 taken 249 times.
✓ Branch 74 → 78 taken 124 times.
✗ Branch 74 → 79 not taken.
1277 call this%insert_polar_extractor_3forms()
210 end select
211
212
4/4
✓ Branch 79 → 80 taken 858 times.
✓ Branch 79 → 85 taken 419 times.
✓ Branch 80 → 81 taken 513 times.
✓ Branch 80 → 85 taken 345 times.
1277 if (this%enriched .and. this%regularity > 0) then
213 ! The backward part of commutativity (weak part of Hodge decomposition)
214 118 select case (space%m)
215 case (1)
216 118 call this%insert_polar_extractor_1forms_enrich()
217 case (2)
218 106 call this%insert_polar_extractor_2forms_enrich()
219 case (3)
220
4/4
✓ Branch 81 → 82 taken 118 times.
✓ Branch 81 → 83 taken 106 times.
✓ Branch 81 → 84 taken 55 times.
✓ Branch 81 → 85 taken 234 times.
513 call this%insert_polar_extractor_3forms_enrich()
221 end select
222 end if
223
224
1/2
✓ Branch 85 → 86 taken 1277 times.
✗ Branch 85 → 106 not taken.
1277 if (this%orthogonal) then
225 ! The 3D polar forms are Kronecker products of identity matrix (z) with 2D polar forms, so 2D orthogonalization is sufficient
226 ! NOTE: E1 and E2 are assumed to act on diffent components of the m-form, hence orthogonalization is independent
227
1/2
✓ Branch 86 → 87 taken 1277 times.
✗ Branch 86 → 96 not taken.
1277 if (associated(this%E1)) then
228 1277 call this%shmem_E1%fence()
229
2/4
✓ Branch 88 → 89 taken 1277 times.
✗ Branch 88 → 94 not taken.
✗ Branch 91 → 92 not taken.
✓ Branch 91 → 94 taken 1277 times.
1277 if (this%shmem_E1%leader()) call orthogonalize(this%E1)
230 1277 call this%shmem_E1%fence()
231 end if
232
2/2
✓ Branch 96 → 97 taken 379 times.
✓ Branch 96 → 106 taken 898 times.
1277 if (associated(this%E2)) then
233 379 call this%shmem_E2%fence()
234
2/4
✓ Branch 98 → 99 taken 379 times.
✗ Branch 98 → 104 not taken.
✗ Branch 101 → 102 not taken.
✓ Branch 101 → 104 taken 379 times.
379 if (this%shmem_E2%leader()) call orthogonalize(this%E2)
235 379 call this%shmem_E2%fence()
236 end if
237 end if
238
239 533 select case (space%m)
240 case (0)
241 533 call this%finalize_polar_extractor_0forms()
242 case (1)
243 371 call this%finalize_polar_extractor_1forms()
244 case (2)
245 249 call this%finalize_polar_extractor_2forms()
246 case (3)
247
4/5
✓ Branch 106 → 107 taken 533 times.
✓ Branch 106 → 108 taken 371 times.
✓ Branch 106 → 109 taken 249 times.
✓ Branch 106 → 110 taken 124 times.
✗ Branch 106 → 111 not taken.
1277 call this%finalize_polar_extractor_3forms()
248 end select
249
250
1/2
✓ Branch 111 → 112 taken 1277 times.
✗ Branch 111 → 114 not taken.
1277 if (associated(this%E1)) then
251 1277 call this%shmem_E1%destroy()
252 1277 nullify (this%E1)
253 end if
254
2/2
✓ Branch 114 → 115 taken 379 times.
✓ Branch 114 → 117 taken 898 times.
1277 if (associated(this%E2)) then
255 379 call this%shmem_E2%destroy()
256 379 nullify (this%E2)
257 end if
258
259 end subroutine init_mform_polar
260
261 !> @brief Initialize the shared memory windows for the extractor matrices
262 !>
263 !> @param[inout] this The MFormPolarRegularity object to initialize
264 !> @param[inout] E4d The 4D extractor matrix to be initialized as a shared memory window
265 !> @param[inout] E_window The SharedMemoryWindow object managing the shared memory for E4d
266 !> @param[in] nr_polar_splines The number of polar splines used in the extractor
267 1845 subroutine init_shmem_window(this, E4d, E_window, nr_polar_splines)
268 use m_common, only: cartesian_to_linear
269 use m_tensorprod, only: size
270
271 implicit none
272
273 class(MFormPolarRegularity), intent(inout) :: this
274 real(wp), pointer :: E4d(:, :, :, :)
275 type(SharedMemoryWindow), intent(inout) :: E_window
276 integer, intent(in) :: nr_polar_splines
277
278 integer :: comm_color, comm_key
279 integer :: nr_x, nr_y, nr_blocks, nr_elements
280 1845 integer, allocatable :: cart_inds(:), cart_size(:), non_cart_inds(:), non_cart_size(:)
281 integer :: ierr
282 logical :: mask(3)
283
284
2/2
✓ Branch 2 → 3 taken 1277 times.
✓ Branch 2 → 19 taken 568 times.
1845 if (this%comm_window == -1) then
285 ! Create a shared memory sub communicator for the ranks involved in the constraint that share mpi unknowns
286 1277 mask(1) = this%space%tp_spaces(1)%domain%is_shared_memory(1)
287 1277 mask(2) = this%space%tp_spaces(1)%domain%is_shared_memory(2)
288 1277 mask(3) = .false.
289
290 1277 cart_inds = pack(this%space%tp_spaces(1)%domain%my_subinterval_ijk, mask=mask)
291 1277 cart_size = pack(this%space%tp_spaces(1)%domain%nr_subintervals, mask=mask)
292
2/2
✓ Branch 6 → 7 taken 3831 times.
✓ Branch 6 → 8 taken 1277 times.
5108 non_cart_inds = pack(this%space%tp_spaces(1)%domain%my_subinterval_ijk, mask=.not. mask)
293
2/2
✓ Branch 10 → 11 taken 3831 times.
✓ Branch 10 → 12 taken 1277 times.
5108 non_cart_size = pack(this%space%tp_spaces(1)%domain%nr_subintervals, mask=.not. mask)
294
295 1277 comm_color = cartesian_to_linear(non_cart_inds, non_cart_size) ! Same color ends up on same comm
296 1277 comm_key = cartesian_to_linear(cart_inds, cart_size) ! Key to order ranks within comm
297
298
1/2
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 1277 times.
1277 PetscCallMPI(MPI_Comm_split(this%space%tp_spaces(1)%domain%comm, comm_color, comm_key, this%comm_window, ierr))
299
1/2
✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 1277 times.
1277 PetscCallMPI(MPI_Comm_rank(this%comm_window, this%my_rank_window, ierr))
300 end if
301
302 nr_x = merge(this%regularity + 1, 0, this%space%tp_spaces(1)%domain%is_shared_memory(1) .or. &
303
2/4
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 21 taken 1845 times.
✓ Branch 21 → 20 taken 1845 times.
✗ Branch 21 → 22 not taken.
1845 this%space%tp_spaces(1)%domain%my_subinterval_ijk(1) == 0)
304 1845 nr_y = size(this%space%tp_spaces(1), 2) ! All spaces have same size in y-direction
305 1845 nr_blocks = size(this%space%tp_spaces)
306
307 1845 nr_elements = nr_x * nr_y * nr_blocks * nr_polar_splines
308
309 1845 call E_window%init(nr_elements, this%comm_window, this%my_rank_window)
310 1845 E4d(0:nr_x - 1, 0:nr_y - 1, 1:nr_blocks, 1:nr_polar_splines) => E_window%window(:)
311
9/10
✓ Branch 23 → 24 taken 1845 times.
✗ Branch 23 → 35 not taken.
✓ Branch 24 → 25 taken 7529 times.
✓ Branch 24 → 35 taken 1845 times.
✓ Branch 26 → 27 taken 14509 times.
✓ Branch 26 → 34 taken 7529 times.
✓ Branch 28 → 29 taken 425714 times.
✓ Branch 28 → 33 taken 14509 times.
✓ Branch 30 → 31 taken 1964634 times.
✓ Branch 30 → 32 taken 425714 times.
2414231 if (E_window%leader()) E4d = 0._wp
312 1845 call E_window%fence()
313
314
8/16
✓ Branch 36 → 37 taken 1277 times.
✓ Branch 36 → 38 taken 568 times.
✓ Branch 38 → 39 taken 1277 times.
✓ Branch 38 → 40 taken 568 times.
✓ Branch 40 → 41 taken 1277 times.
✓ Branch 40 → 42 taken 568 times.
✓ Branch 42 → 43 taken 1277 times.
✓ Branch 42 → 59 taken 568 times.
✗ Branch 44 → 45 not taken.
✗ Branch 44 → 46 not taken.
✗ Branch 48 → 49 not taken.
✗ Branch 48 → 50 not taken.
✗ Branch 52 → 53 not taken.
✗ Branch 52 → 54 not taken.
✗ Branch 56 → 57 not taken.
✗ Branch 56 → 60 not taken.
1845 end subroutine init_shmem_window
315
316 !> @brief Apply the polar regularity constraint to an m-form (overload of MFormConstraintLocal::apply)
317 !>
318 !> @param[in] this The constraint to apply
319 !> @param[out] vout The MForm to which the constraint has been applied
320 !> @param[in] vin The MForm to which the constraint is applied
321
4/4
✓ Branch 2 → 3 taken 388 times.
✓ Branch 2 → 4 taken 323 times.
✓ Branch 4 → 5 taken 388 times.
✓ Branch 4 → 6 taken 323 times.
711 subroutine apply_mform_polar(this, vout, vin)
322 use m_mform_basis, only: MFormFun
323 implicit none
324
325 class(MFormPolarRegularity), intent(in) :: this
326 type(MFormFun), intent(out) :: vout
327 type(MFormFun), intent(in) :: vin
328
329 integer :: blk, local_idx, local_row, local_col
330 integer :: i, j, k, row_start, row_end
331 real(wp) :: local_val
332 real(wp), allocatable :: reduced_coeffs(:), my_reduced_coeffs(:)
333 integer :: nr_cols, ierr
334
335 ! The columns of the extraction operator are orthonormal, so we can use the following projection
336 !> vout = E * (E^T * vin)
337 !>
338 !> This is done in two steps:
339 !> 1. Compute reduced_coeffs = E^T * vin (project vin onto the constrained subspace)
340 !> 2. Compute vout = E * reduced_coeffs (reconstruct in the full space)
341
342 ! Make a deep copy
343 711 vout = vin
344
345
1/2
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 711 times.
711 if (.not. allocated(this%extractor)) error stop "MFormPolarRegularity::apply_mform_polar: extractor not allocated"
346
1/2
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 711 times.
711 if (.not. this%orthogonal) error stop "MFormPolarRegularity::apply_mform_polar: only orthogonal extraction can be used in apply"
347
1/2
✗ Branch 11 → 12 not taken.
✓ Branch 11 → 13 taken 711 times.
711 if (this%comm_window == -1) error stop "MFormPolarRegularity::apply_mform_polar: shared memory window communicator not"// &
348 " initialized"
349
350 ! Apply projection using extraction operator: vout = E * (E^T * vin)
351 ! Only apply on the rank at the polar singularity (my_rank_i == 0)
352 ! TODO node-level-mpi: requires shared memory window with sparse extractor matrices
353 ! Set to zero where the constraint acts
354
2/2
✓ Branch 14 → 15 taken 1797 times.
✓ Branch 14 → 20 taken 711 times.
2508 do blk = 1, size(vout%tp_funs)
355
2/2
✓ Branch 16 → 17 taken 145052 times.
✓ Branch 16 → 18 taken 1797 times.
147560 do local_idx = 0, size(this%index_list(blk)) - 1
356 145052 i = this%index_list(blk)%i(local_idx)
357 145052 j = this%index_list(blk)%j(local_idx)
358 145052 k = this%index_list(blk)%k(local_idx)
359
360 146849 vout%tp_funs(blk)%data(i, j, k) = 0._wp
361 end do
362 end do
363
364 711 nr_cols = this%extractor(1)%nr_cols
365
8/12
✓ Branch 21 → 22 taken 439 times.
✓ Branch 21 → 23 taken 272 times.
✓ Branch 23 → 22 taken 272 times.
✗ Branch 23 → 24 not taken.
✓ Branch 24 → 25 taken 711 times.
✗ Branch 24 → 26 not taken.
✓ Branch 26 → 27 taken 272 times.
✓ Branch 26 → 28 taken 439 times.
✗ Branch 28 → 29 not taken.
✓ Branch 28 → 30 taken 711 times.
✗ Branch 30 → 31 not taken.
✓ Branch 30 → 32 taken 711 times.
2133 allocate (my_reduced_coeffs(0:nr_cols - 1))
366
2/2
✓ Branch 32 → 33 taken 4884 times.
✓ Branch 32 → 34 taken 711 times.
5595 my_reduced_coeffs = 0._wp
367
368 ! Step 1: Compute E^T * vin (transpose multiply)
369 ! Loop over all input blocks that contribute
370
2/2
✓ Branch 35 → 36 taken 1797 times.
✓ Branch 35 → 44 taken 711 times.
2508 do blk = 1, size(vin%tp_funs)
371
2/2
✓ Branch 37 → 38 taken 145052 times.
✓ Branch 37 → 42 taken 1797 times.
147560 do local_row = 0, this%extractor(blk)%nr_rows - 1
372 145052 row_start = this%extractor(blk)%row_starts_at_nz(local_row)
373 145052 row_end = this%extractor(blk)%row_starts_at_nz(local_row + 1) - 1
374
375 145052 i = this%index_list(blk)%i(local_row)
376 145052 j = this%index_list(blk)%j(local_row)
377 145052 k = this%index_list(blk)%k(local_row)
378
379
2/2
✓ Branch 39 → 40 taken 145052 times.
✓ Branch 39 → 41 taken 453180 times.
600029 do local_idx = row_start, row_end
380 453180 local_col = this%extractor(blk)%col_idx(local_idx)
381 453180 local_val = this%extractor(blk)%values(local_idx)
382
383 ! E^T multiply: reduced_coeffs(col) += E(col, row) * vin(row)
384 598232 my_reduced_coeffs(local_col) = my_reduced_coeffs(local_col) + local_val * vin%tp_funs(blk)%data(i, j, k)
385 end do
386 end do
387 end do
388
389
5/8
✓ Branch 45 → 46 taken 439 times.
✓ Branch 45 → 47 taken 272 times.
✓ Branch 47 → 46 taken 272 times.
✗ Branch 47 → 48 not taken.
✗ Branch 48 → 49 not taken.
✓ Branch 48 → 50 taken 711 times.
✗ Branch 50 → 51 not taken.
✓ Branch 50 → 52 taken 711 times.
1422 allocate (reduced_coeffs(0:nr_cols - 1))
390 ! TODO shared memory equivalent? yes, use a fenced shmem window
391
1/2
✗ Branch 53 → 54 not taken.
✓ Branch 53 → 57 taken 711 times.
711 PetscCallMPI(MPI_Allreduce(my_reduced_coeffs, reduced_coeffs, nr_cols, MPI_WP, MPI_SUM, this%comm_window, ierr))
392
393 ! Step 2: Compute E * reduced_coeffs (forward multiply)
394
2/2
✓ Branch 58 → 59 taken 1797 times.
✓ Branch 58 → 67 taken 711 times.
2508 do blk = 1, size(vout%tp_funs)
395
2/2
✓ Branch 60 → 61 taken 145052 times.
✓ Branch 60 → 65 taken 1797 times.
147560 do local_row = 0, this%extractor(blk)%nr_rows - 1
396 145052 row_start = this%extractor(blk)%row_starts_at_nz(local_row)
397 145052 row_end = this%extractor(blk)%row_starts_at_nz(local_row + 1) - 1
398
399 145052 i = this%index_list(blk)%i(local_row)
400 145052 j = this%index_list(blk)%j(local_row)
401 145052 k = this%index_list(blk)%k(local_row)
402
403
2/2
✓ Branch 62 → 63 taken 145052 times.
✓ Branch 62 → 64 taken 453180 times.
600029 do local_idx = row_start, row_end
404 453180 local_col = this%extractor(blk)%col_idx(local_idx)
405 453180 local_val = this%extractor(blk)%values(local_idx)
406
407 ! E multiply: vout(row) += E(row, col) * reduced_coeffs(col)
408 vout%tp_funs(blk)%data(i, j, k) = &
409 598232 vout%tp_funs(blk)%data(i, j, k) + local_val * reduced_coeffs(local_col)
410 end do
411 end do
412 end do
413
414
1/2
✗ Branch 68 → 69 not taken.
✓ Branch 68 → 70 taken 711 times.
711 deallocate (reduced_coeffs)
415
1/2
✗ Branch 70 → 71 not taken.
✓ Branch 70 → 72 taken 711 times.
711 deallocate (my_reduced_coeffs)
416
417 711 call vout%distribute()
418 end subroutine apply_mform_polar
419
420 !> @brief Initialize the 0-form extractor for the polar regularity constraint
421 !>
422 !> @param[inout] this The MFormPolarRegularity object to initialize
423 533 subroutine insert_polar_extractor_0forms(this)
424 use m_tensorprod, only: size
425 use m_tensorprod_indices, only: TensorProdIndices
426 implicit none
427
428 class(MFormPolarRegularity), intent(inout) :: this
429
430 type(TensorProdIndices) :: inds
431
432 integer :: i, j, i_min, i_max, j_min, j_max, k_min, k_max
433 integer :: gamma, alpha, polar_index, nr_polar_splines, gamma_min, gamma_max
434 integer :: s1
435 real(wp) :: coeff
436
437 533 gamma_min = 0
438 533 gamma_max = this%regularity
439 533 nr_polar_splines = get_nr_polar_zeroforms(gamma_max)
440
441 !> Initialize the indices on which this extractor acts
442 533 i_min = max(0, this%space%tp_spaces(1)%rank_resp_bspline%i0)
443 533 i_max = min(this%regularity, this%space%tp_spaces(1)%rank_resp_bspline%i1)
444 533 j_min = this%space%tp_spaces(1)%rank_resp_bspline%j0; j_max = this%space%tp_spaces(1)%rank_resp_bspline%j1
445 533 k_min = this%space%tp_spaces(1)%rank_resp_bspline%k0; k_max = this%space%tp_spaces(1)%rank_resp_bspline%k1
446 533 call inds%init(i_min, i_max, j_min, j_max, k_min, k_max, this%space%tp_spaces(1)%spaces)
447 533 call this%index_list(1)%init(inds, this%space%tp_spaces(1)%rank_resp_bspline, this%space%tp_spaces(1)%rank_l0)
448
2/2
✓ Branch 5 → 6 taken 2132 times.
✓ Branch 5 → 7 taken 533 times.
2665 this%row0_2d = 0
449
450 !> Initialize the radial and angular B-splines
451 call this%compute_1d_bsplines(this%space%tp_spaces(1)%spaces(1), this%space%tp_spaces(1)%spaces(2), gamma_min, gamma_max, &
452 533 radial_truncate=this%regularity)
453
454 533 call this%init_shmem_window(this%E1, this%shmem_E1, nr_polar_splines)
455 s1 = 1
456
2/2
✓ Branch 9 → 10 taken 15548 times.
✓ Branch 9 → 25 taken 533 times.
16081 do j = j_min, j_max
457
2/2
✓ Branch 11 → 12 taken 44506 times.
✓ Branch 11 → 23 taken 15548 times.
60587 do i = i_min, i_max
458 polar_index = 0
459
2/2
✓ Branch 13 → 14 taken 187898 times.
✓ Branch 13 → 21 taken 44506 times.
247952 do gamma = gamma_min, gamma_max
460
1/2
✓ Branch 14 → 15 taken 187898 times.
✗ Branch 14 → 19 not taken.
232404 do alpha = -gamma, gamma, 2
461
2/2
✓ Branch 16 → 17 taken 187898 times.
✓ Branch 16 → 18 taken 405904 times.
593802 this%E1(i, j, s1, polar_index + 1) = this%radial_bsplines(gamma)%data(i) * this%angular_bsplines(alpha)%data(j)
462 187898 polar_index = polar_index + 1
463 end do
464 end do
465 end do
466 end do
467
468 533 end subroutine insert_polar_extractor_0forms
469
470 !> @brief Finalize the 0-form extractor for the polar regularity constraint
471 !>
472 !> @param[inout] this The MFormPolarRegularity object to initialize
473 533 subroutine finalize_polar_extractor_0forms(this)
474 implicit none
475
476 class(MFormPolarRegularity), intent(inout) :: this
477
478 integer :: nr_nonzeros
479 integer :: i, j, k, row_3d, col, polar_index, s1
480 integer :: i_min, i_max, j_min, j_max, k_min, k_max, nr_polar_splines
481 real(wp) :: coeff
482
483 533 call this%index_list(1)%get_bounds(i_min, i_max, j_min, j_max, k_min, k_max)
484
485 533 nr_polar_splines = size(this%E1, 4)
486
487 !> Insert the tensor product of radial and angular B-splines into the extractor
488 533 nr_nonzeros = size(this%index_list(1)) * nr_polar_splines
489 call this%extractor(1)%init(nr_rows=size(this%index_list(1)), nr_cols=nr_polar_splines * (1 + k_max - k_min), &
490 533 nr_nonzero=nr_nonzeros)
491
492 s1 = 1
493
2/2
✓ Branch 5 → 6 taken 1385 times.
✓ Branch 5 → 22 taken 533 times.
1918 do k = k_min, k_max
494
2/2
✓ Branch 7 → 8 taken 36880 times.
✓ Branch 7 → 20 taken 1385 times.
38798 do j = j_min, j_max
495
2/2
✓ Branch 9 → 10 taken 122670 times.
✓ Branch 9 → 18 taken 36880 times.
160935 do i = i_min, i_max
496 122670 row_3d = this%index_list(1)%get_local_index(i, j, k)
497
498
2/2
✓ Branch 12 → 13 taken 1548734 times.
✓ Branch 12 → 16 taken 122670 times.
1708284 do polar_index = 0, nr_polar_splines - 1
499 1548734 col = polar_index + (k - k_min) * nr_polar_splines
500 1548734 coeff = this%E1(i, j, s1, polar_index + 1)
501
2/2
✓ Branch 13 → 11 taken 366282 times.
✓ Branch 13 → 14 taken 1182452 times.
1671404 if (abs(coeff) > COEFF_THR) call this%extractor(1)%insert(row_3d, col, coeff)
502 end do
503 end do
504 end do
505 end do
506
507 533 call this%extractor(1)%finalize()
508 533 end subroutine finalize_polar_extractor_0forms
509
510 !> @brief Initialize the 1-form extractor for the polar regularity constraint using gradients of polar zero-forms
511 !>
512 !> @param[inout] this The MFormPolarRegularity object to initialize
513 371 subroutine insert_polar_extractor_1forms(this)
514 use m_tensorprod, only: size
515 use m_qr_factorisation, only: orthogonalize
516 use m_tensorprod_indices, only: TensorProdIndices
517 implicit none
518
519 class(MFormPolarRegularity), intent(inout) :: this
520
521 type(TensorProdIndices) :: inds
522
523 integer :: nr_polar_oneforms_xy, nr_polar_splines_z
524 integer :: nr_bspline_y
525 integer :: i, j, k, row, s1, polar_index, col, row0_2d(4), j_plus_one
526 integer :: i_min(3), i_max(3), j_min(3), j_max(3), k_min(3), k_max(3)
527 integer :: gamma, alpha, gamma_min(3), gamma_max
528 real(wp) :: coeff
529
530 ! NOTE: the y,z spaces are periodic and therefore the number of basis functions is the same
531 371 nr_bspline_y = size(this%space%tp_spaces(1), 2)
532
533 371 gamma_min(1) = 1
534 371 gamma_min(2) = 1
535 gamma_min(3) = 0
536 371 gamma_max = this%regularity
537
538 371 nr_polar_oneforms_xy = get_nr_polar_zeroforms(gamma_max) - get_nr_polar_zeroforms(0)
539 742 nr_polar_splines_z = merge(get_nr_polar_zeroforms(gamma_max), 0, this%space%nr_components == 3)
540
541 !> Initialize the indices on which this extractor acts
542 371 row0_2d(1) = 0
543
2/2
✓ Branch 5 → 6 taken 1113 times.
✓ Branch 5 → 11 taken 371 times.
1484 do s1 = 1, 3
544 1113 i_min(s1) = max(0, this%space%tp_spaces(s1)%rank_resp_bspline%i0)
545
2/2
✓ Branch 6 → 7 taken 742 times.
✓ Branch 6 → 8 taken 371 times.
1113 i_max(s1) = this%regularity - merge(1, 0, s1 == 1)
546 1113 i_max(s1) = min(i_max(s1), this%space%tp_spaces(s1)%rank_resp_bspline%i1)
547 1113 j_min(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%j0
548 1113 j_max(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%j1
549 1113 k_min(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%k0
550 1113 k_max(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%k1
551 1113 call inds%init(i_min(s1), i_max(s1), j_min(s1), j_max(s1), k_min(s1), k_max(s1), this%space%tp_spaces(s1)%spaces)
552 1113 call this%index_list(s1)%init(inds, this%space%tp_spaces(s1)%rank_resp_bspline, this%space%tp_spaces(s1)%rank_l0)
553 1484 row0_2d(s1 + 1) = row0_2d(s1) + size_2d(this%index_list(s1))
554 end do
555
2/2
✓ Branch 13 → 14 taken 1484 times.
✓ Branch 13 → 15 taken 371 times.
1855 this%row0_2d = row0_2d
556
557 !> Initialize the radial and angular B-splines (data members of this object, using space0!)
558 call this%compute_1d_bsplines(this%space0%tp_spaces(1)%spaces(1), this%space0%tp_spaces(1)%spaces(2), 0, this%regularity, &
559 371 radial_truncate=this%regularity)
560
561 ! xy-components
562 371 call this%init_shmem_window(this%E1, this%shmem_E1, nr_polar_oneforms_xy)
563
2/2
✓ Branch 17 → 18 taken 742 times.
✓ Branch 17 → 38 taken 371 times.
1113 do s1 = 1, 2
564
2/2
✓ Branch 19 → 20 taken 21960 times.
✓ Branch 19 → 36 taken 742 times.
23073 do j = j_min(s1), j_max(s1)
565
2/2
✓ Branch 21 → 22 taken 38664 times.
✓ Branch 21 → 34 taken 21960 times.
61366 do i = i_min(s1), i_max(s1)
566 polar_index = 0
567 ! NOTE: we skip gamma == 0 because it's gradient is zero
568
2/2
✓ Branch 23 → 24 taken 99110 times.
✓ Branch 23 → 32 taken 38664 times.
159734 do gamma = gamma_min(s1), gamma_max
569
1/2
✗ Branch 24 → 25 not taken.
✓ Branch 24 → 26 taken 99110 times.
137774 do alpha = -gamma, gamma, 2
570 142872 select case (s1)
571 case (1)
572 coeff = (this%radial_bsplines(gamma)%data(i + 1) - this%radial_bsplines(gamma)%data(i)) * &
573 142872 this%angular_bsplines(alpha)%data(j)
574 case (2)
575 184952 j_plus_one = mod(j + 1, nr_bspline_y)
576 coeff = this%radial_bsplines(gamma)%data(i) * &
577
2/2
✓ Branch 26 → 27 taken 142872 times.
✓ Branch 26 → 28 taken 184952 times.
327824 (this%angular_bsplines(alpha)%data(j_plus_one) - this%angular_bsplines(alpha)%data(j))
578 end select
579
2/2
✓ Branch 29 → 26 taken 228714 times.
✓ Branch 29 → 30 taken 99110 times.
327824 this%E1(i, j, s1, polar_index + 1) = coeff
580 99110 polar_index = polar_index + 1
581 end do
582 end do
583 end do
584 end do
585 end do
586
587 ! z-component
588
2/2
✓ Branch 39 → 40 taken 130 times.
✓ Branch 39 → 56 taken 241 times.
371 if (this%space%nr_components == 3) then
589 130 call this%init_shmem_window(this%E2, this%shmem_E2, nr_polar_splines_z)
590
591 s1 = 3
592
2/2
✓ Branch 42 → 43 taken 3772 times.
✓ Branch 42 → 55 taken 130 times.
3902 do j = j_min(s1), j_max(s1)
593
2/2
✓ Branch 44 → 45 taken 10528 times.
✓ Branch 44 → 53 taken 3772 times.
14430 do i = i_min(s1), i_max(s1)
594 polar_index = 0
595
2/2
✓ Branch 45 → 46 taken 37956 times.
✓ Branch 45 → 51 taken 10528 times.
52256 do gamma = gamma_min(s1), gamma_max
596
1/2
✗ Branch 46 → 47 not taken.
✓ Branch 46 → 48 taken 37956 times.
48484 do alpha = -gamma, gamma, 2
597 this%E2(i, j, s1, polar_index + 1) = &
598
2/2
✓ Branch 48 → 48 taken 62060 times.
✓ Branch 48 → 49 taken 37956 times.
100016 this%radial_bsplines(gamma)%data(i) * this%angular_bsplines(alpha)%data(j)
599 37956 polar_index = polar_index + 1
600 end do
601 end do
602 end do
603 end do
604 end if
605
606 371 end subroutine insert_polar_extractor_1forms
607
608 !> @brief Enrich 1-form basis for the polar regularity constraint using weak curls of polar two-forms
609 !>
610 !> @param[inout] this The MFormPolarEnriched object to initialize
611 118 subroutine insert_polar_extractor_1forms_enrich(this)
612 use m_bspline_basis, only: BSplineFun, BSplineSpace, evaluate
613 use m_bspline_linalg, only: l2_projection, biorthogonalize
614 use m_common, only: zero_function, wp, pi, VERBOSITY_WARN_ALWAYS
615 use m_coord_transform, only: CylinderTransform
616 use m_mform_basis, only: MFormSpace, get_zero_form_space, evaluate, MFormFun
617 use m_mform_constraint_boundary, only: MFormDirichlet, MFormNeumann
618 use m_mform_derham, only: DeRhamSequence
619 use m_mform_linalg, only: inner_product, l2_norm, l2_error, solve_curlA_equals_B
620 use m_mform_solver, only: GenericSolver
621 use m_mform_matrix, only: DiffDiffMatrix, MassMatrix
622 use m_sparsemat, only: SparseMat
623 use m_tensorprod, only: size
624 use m_domain_decomp, only: DomainDecomp, memory_layout_convert
625 implicit none
626
627 class(MFormPolarRegularity), intent(inout) :: this
628
629 type(CylinderTransform) :: disk
630 118 type(GenericSolver) :: divgrad_solver_2d
631 118 type(DiffDiffMatrix) :: divgrad_matrix_2d
632 type(DomainDecomp) :: domain_2d
633
6/6
✓ Branch 3 → 4 taken 354 times.
✓ Branch 3 → 5 taken 118 times.
✓ Branch 280 → 281 taken 354 times.
✓ Branch 280 → 284 taken 118 times.
✓ Branch 281 → 282 taken 240 times.
✓ Branch 281 → 283 taken 114 times.
944 type(MFormSpace) :: derham_2d(0:2)
634 118 type(MFormFun) :: f, tmp
635 type(MFormFun), allocatable :: polar_oneforms(:)
636 118 type(MFormPolarRegularity) :: bounded_oneform_constraint
637
638 type(BSplineSpace) :: bspline_x, bspline_y
639
640 118 real(wp), pointer :: E1_curl(:, :, :, :)
641 type(SharedMemoryWindow) :: shmem_E1_curl
642
643 integer :: comm_2d, ierr, comm_color, comm_key
644 integer :: nr_enriched_basis, nr_gradient_basis, gamma, alpha, gamma_min, gamma_max
645 integer :: s1, nr_nonzeros(3), i, j, row, polar_index, nr_x, nr_bspline_y
646 integer :: i_min(3), i_max(3), j_min(3), j_max(3)
647
648 gamma_min = 0
649
650 ! The pushforward yields extra factor xp, and the 3rd component of the 2-form space has degree reduced by 1
651 118 gamma_max = this%regularity - 2
652 118 if (gamma_max < gamma_min) return
653
654 80 disk = CylinderTransform(radius=1.0_wp)
655
656 ! Create a 2D domain decomposition for the disk cross-section, respecting the original 3D decomposition
657 80 comm_color = this%space%tp_spaces(1)%domain%my_subinterval_ijk(3)
658 comm_key = this%space%tp_spaces(1)%domain%my_subinterval_ijk(1) + this%space%tp_spaces(1)%domain%my_subinterval_ijk(2) * &
659 80 this%space%tp_spaces(1)%domain%nr_subintervals(1)
660
1/2
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 80 times.
80 PetscCallMPI(MPI_Comm_split(this%space%tp_spaces(1)%domain%comm, comm_color, comm_key, comm_2d, ierr))
661 domain_2d = DomainDecomp(&
662 & Nx=this%space%tp_spaces(1)%domain%nr_subintervals(1), &
663 & Ny=this%space%tp_spaces(1)%domain%nr_subintervals(2), &
664 & memory_layout_x=memory_layout_convert(this%space%tp_spaces(1)%domain%memory_layout(1)), &
665 & memory_layout_y=memory_layout_convert(this%space%tp_spaces(1)%domain%memory_layout(2)), &
666 80 comm=comm_2d, flat_mpi=this%space%tp_spaces(1)%domain%flat_mpi, dimensionality=2)
667
668 ! Get the 2D DeRham sequence corresponding to the disk cross-section
669
3/4
✓ Branch 15 → 16 taken 240 times.
✓ Branch 15 → 19 taken 80 times.
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 240 times.
320 derham_2d = DeRhamSequence(this%space0%tp_spaces(1)%spaces(1), this%space0%tp_spaces(1)%spaces(2), domain=domain_2d)
670
671 ! Set up the curl-curl and divgrad solver in 2D
672 80 call divgrad_matrix_2d%init(derham_2d(0), disk)
673 call divgrad_solver_2d%init(divgrad_matrix_2d, constraint1=MFormNeumann(), constraint2=MFormPolarBoundedness(), &
674
5/36
✓ Branch 21 → 22 taken 80 times.
✗ Branch 21 → 23 not taken.
✗ Branch 24 → 25 not taken.
✓ Branch 24 → 26 taken 80 times.
✗ Branch 26 → 27 not taken.
✓ Branch 26 → 39 taken 80 times.
✗ Branch 28 → 29 not taken.
✗ Branch 28 → 38 not taken.
✗ Branch 29 → 30 not taken.
✗ Branch 29 → 31 not taken.
✗ Branch 31 → 32 not taken.
✗ Branch 31 → 33 not taken.
✗ Branch 33 → 34 not taken.
✗ Branch 33 → 35 not taken.
✗ Branch 35 → 36 not taken.
✗ Branch 35 → 37 not taken.
✗ Branch 39 → 40 not taken.
✓ Branch 39 → 50 taken 80 times.
✗ Branch 41 → 42 not taken.
✗ Branch 41 → 49 not taken.
✗ Branch 42 → 43 not taken.
✗ Branch 42 → 44 not taken.
✗ Branch 44 → 45 not taken.
✗ Branch 44 → 46 not taken.
✗ Branch 46 → 47 not taken.
✗ Branch 46 → 48 not taken.
✗ Branch 50 → 51 not taken.
✓ Branch 50 → 61 taken 80 times.
✗ Branch 52 → 53 not taken.
✗ Branch 52 → 60 not taken.
✗ Branch 53 → 54 not taken.
✗ Branch 53 → 55 not taken.
✗ Branch 55 → 56 not taken.
✗ Branch 55 → 57 not taken.
✗ Branch 57 → 58 not taken.
✗ Branch 57 → 59 not taken.
160 coord_transform=disk)
675
676
7/48
✗ Branch 61 → 62 not taken.
✓ Branch 61 → 63 taken 80 times.
✗ Branch 63 → 64 not taken.
✓ Branch 63 → 76 taken 80 times.
✗ Branch 65 → 66 not taken.
✗ Branch 65 → 75 not taken.
✗ Branch 66 → 67 not taken.
✗ Branch 66 → 68 not taken.
✗ Branch 68 → 69 not taken.
✗ Branch 68 → 70 not taken.
✗ Branch 70 → 71 not taken.
✗ Branch 70 → 72 not taken.
✗ Branch 72 → 73 not taken.
✗ Branch 72 → 74 not taken.
✗ Branch 76 → 77 not taken.
✓ Branch 76 → 87 taken 80 times.
✗ Branch 78 → 79 not taken.
✗ Branch 78 → 86 not taken.
✗ Branch 79 → 80 not taken.
✗ Branch 79 → 81 not taken.
✗ Branch 81 → 82 not taken.
✗ Branch 81 → 83 not taken.
✗ Branch 83 → 84 not taken.
✗ Branch 83 → 85 not taken.
✗ Branch 87 → 88 not taken.
✓ Branch 87 → 98 taken 80 times.
✗ Branch 89 → 90 not taken.
✗ Branch 89 → 97 not taken.
✗ Branch 90 → 91 not taken.
✗ Branch 90 → 92 not taken.
✗ Branch 92 → 93 not taken.
✗ Branch 92 → 94 not taken.
✗ Branch 94 → 95 not taken.
✗ Branch 94 → 96 not taken.
✗ Branch 98 → 99 not taken.
✓ Branch 98 → 105 taken 80 times.
✗ Branch 100 → 101 not taken.
✗ Branch 100 → 104 not taken.
✗ Branch 101 → 102 not taken.
✗ Branch 101 → 103 not taken.
✗ Branch 105 → 106 not taken.
✓ Branch 105 → 112 taken 80 times.
✗ Branch 107 → 108 not taken.
✗ Branch 107 → 111 not taken.
✗ Branch 108 → 109 not taken.
✗ Branch 108 → 110 not taken.
✗ Branch 112 → 113 not taken.
✓ Branch 112 → 114 taken 80 times.
80 bounded_oneform_constraint = MFormPolarBoundedness()
677 80 call bounded_oneform_constraint%init(derham_2d(1), coord_transform=disk)
678
679 80 nr_gradient_basis = size(this%E1, 4)
680 80 nr_enriched_basis = get_nr_polar_zeroforms(gamma_max) - get_nr_polar_zeroforms(gamma_min - 1)
681
6/10
✓ Branch 117 → 116 taken 80 times.
✗ Branch 117 → 118 not taken.
✓ Branch 118 → 119 taken 80 times.
✗ Branch 118 → 120 not taken.
✗ Branch 120 → 121 not taken.
✓ Branch 120 → 122 taken 80 times.
✗ Branch 122 → 123 not taken.
✓ Branch 122 → 124 taken 80 times.
✓ Branch 125 → 126 taken 284 times.
✓ Branch 125 → 127 taken 80 times.
524 allocate (polar_oneforms(0:nr_enriched_basis - 1))
682
683
2/2
✓ Branch 128 → 129 taken 240 times.
✓ Branch 128 → 131 taken 80 times.
320 do s1 = 1, 3
684 320 call this%index_list(s1)%get_bounds(i_min(s1), i_max(s1), j_min(s1), j_max(s1))
685 end do
686
687 ! Compute L2 approximation of the radial and angular functions which make up the Zernike basis
688 ! (see also MFormPolarRegularity::compute_1d_bsplines)
689 80 bspline_x = derham_2d(2)%tp_spaces(3)%spaces(1) ! Radial B-spline space for 2-forms
690 80 bspline_y = derham_2d(2)%tp_spaces(3)%spaces(2) ! Angular B-spline space for 2-forms
691
692 ! NB we let gamma_max = gamma_max + 1 here because the radial part of the 2-form has an extra factor r
693 80 call this%compute_1d_bsplines(bspline_x, bspline_y, gamma_min=1, gamma_max=gamma_max + 1, radial_truncate=gamma_max + 1)
694
695 polar_index = 0
696
2/2
✓ Branch 134 → 135 taken 160 times.
✓ Branch 134 → 142 taken 80 times.
240 do gamma = gamma_min, gamma_max
697
1/2
✗ Branch 135 → 136 not taken.
✓ Branch 135 → 137 taken 160 times.
240 do alpha = -gamma, gamma, 2
698 call solve_curlA_equals_B(tmp, this%radial_bsplines(gamma + 1), this%angular_bsplines(alpha), &
699 284 twoform_space=derham_2d(2), cylinder_transform=disk, zeroform_solver=divgrad_solver_2d)
700
701 284 call bounded_oneform_constraint%apply(polar_oneforms(polar_index), tmp)
702
2/2
✓ Branch 139 → 137 taken 124 times.
✓ Branch 139 → 140 taken 160 times.
284 polar_index = polar_index + 1
703 end do
704 end do
705
706 ! Set up the enriched extractor matrices; only for the ranks that are at the polar axis
707 80 call this%init_shmem_window(E1_curl, shmem_E1_curl, nr_enriched_basis)
708
2/2
✓ Branch 144 → 145 taken 160 times.
✓ Branch 144 → 157 taken 80 times.
240 do s1 = 1, 2
709
2/2
✓ Branch 146 → 147 taken 4404 times.
✓ Branch 146 → 155 taken 160 times.
4644 do j = j_min(s1), j_max(s1)
710
2/2
✓ Branch 148 → 149 taken 15518 times.
✓ Branch 148 → 153 taken 4404 times.
20082 do i = i_min(s1), i_max(s1)
711
2/2
✓ Branch 149 → 150 taken 71492 times.
✓ Branch 149 → 151 taken 15518 times.
91414 do polar_index = 0, nr_enriched_basis - 1
712 87010 E1_curl(i, j, s1, polar_index + 1) = polar_oneforms(polar_index)%tp_funs(s1)%data(i, j, 0)
713 end do
714 end do
715 end do
716 end do
717
718 ! Concatenate the shared memory windows
719 80 call this%shmem_E1%cat(shmem_E1_curl)
720 80 nr_bspline_y = size(this%space%tp_spaces(1), 2)
721 nr_x = merge(this%regularity + 1, 0, this%space%tp_spaces(1)%domain%is_shared_memory(1) .or. &
722
2/4
✗ Branch 159 → 160 not taken.
✓ Branch 159 → 161 taken 80 times.
✓ Branch 161 → 160 taken 80 times.
✗ Branch 161 → 162 not taken.
80 this%space%tp_spaces(1)%domain%my_subinterval_ijk(1) == 0)
723 80 this%E1(0:nr_x - 1, 0:nr_bspline_y - 1, 1:3, 1:(nr_gradient_basis + nr_enriched_basis)) => this%shmem_E1%window(:)
724
725 ! Clean up
726 80 call shmem_E1_curl%destroy()
727 80 nullify (E1_curl)
728
729
1/2
✗ Branch 164 → 165 not taken.
✓ Branch 164 → 166 taken 80 times.
80 PetscCallMPI(MPI_Comm_free(comm_2d, ierr))
730 80 call divgrad_matrix_2d%destroy()
731 80 call divgrad_solver_2d%destroy()
732
733 polar_index = 0
734
2/2
✓ Branch 169 → 170 taken 160 times.
✓ Branch 169 → 176 taken 80 times.
240 do gamma = gamma_min, gamma_max
735
1/2
✗ Branch 170 → 171 not taken.
✓ Branch 170 → 172 taken 160 times.
240 do alpha = -gamma, gamma, 2
736 284 call polar_oneforms(polar_index)%destroy()
737
2/2
✓ Branch 173 → 172 taken 124 times.
✓ Branch 173 → 174 taken 160 times.
284 polar_index = polar_index + 1
738 end do
739 end do
740
741
6/10
✓ Branch 177 → 178 taken 80 times.
✗ Branch 177 → 184 not taken.
✓ Branch 178 → 179 taken 284 times.
✓ Branch 178 → 184 taken 80 times.
✓ Branch 179 → 180 taken 284 times.
✗ Branch 179 → 181 not taken.
✗ Branch 181 → 182 not taken.
✓ Branch 181 → 183 taken 284 times.
✗ Branch 184 → 185 not taken.
✓ Branch 184 → 186 taken 80 times.
364 deallocate (polar_oneforms)
742
46/138
✓ Branch 5 → 6 taken 38 times.
✓ Branch 5 → 7 taken 80 times.
✓ Branch 187 → 188 taken 80 times.
✓ Branch 187 → 189 taken 38 times.
✓ Branch 189 → 190 taken 80 times.
✓ Branch 189 → 191 taken 38 times.
✗ Branch 191 → 192 not taken.
✓ Branch 191 → 199 taken 118 times.
✗ Branch 192 → 193 not taken.
✗ Branch 192 → 198 not taken.
✗ Branch 193 → 194 not taken.
✗ Branch 193 → 195 not taken.
✗ Branch 195 → 196 not taken.
✗ Branch 195 → 197 not taken.
✓ Branch 199 → 200 taken 80 times.
✓ Branch 199 → 201 taken 38 times.
✓ Branch 201 → 202 taken 80 times.
✓ Branch 201 → 203 taken 38 times.
✗ Branch 203 → 204 not taken.
✓ Branch 203 → 275 taken 118 times.
✗ Branch 204 → 205 not taken.
✗ Branch 204 → 206 not taken.
✗ Branch 207 → 208 not taken.
✗ Branch 207 → 251 not taken.
✗ Branch 208 → 209 not taken.
✗ Branch 208 → 211 not taken.
✗ Branch 209 → 210 not taken.
✗ Branch 209 → 211 not taken.
✗ Branch 211 → 212 not taken.
✗ Branch 211 → 250 not taken.
✗ Branch 212 → 213 not taken.
✗ Branch 212 → 214 not taken.
✗ Branch 214 → 215 not taken.
✗ Branch 214 → 227 not taken.
✗ Branch 216 → 217 not taken.
✗ Branch 216 → 226 not taken.
✗ Branch 217 → 218 not taken.
✗ Branch 217 → 219 not taken.
✗ Branch 219 → 220 not taken.
✗ Branch 219 → 221 not taken.
✗ Branch 221 → 222 not taken.
✗ Branch 221 → 223 not taken.
✗ Branch 223 → 224 not taken.
✗ Branch 223 → 225 not taken.
✗ Branch 227 → 228 not taken.
✗ Branch 227 → 238 not taken.
✗ Branch 229 → 230 not taken.
✗ Branch 229 → 237 not taken.
✗ Branch 230 → 231 not taken.
✗ Branch 230 → 232 not taken.
✗ Branch 232 → 233 not taken.
✗ Branch 232 → 234 not taken.
✗ Branch 234 → 235 not taken.
✗ Branch 234 → 236 not taken.
✗ Branch 238 → 239 not taken.
✗ Branch 238 → 249 not taken.
✗ Branch 240 → 241 not taken.
✗ Branch 240 → 248 not taken.
✗ Branch 241 → 242 not taken.
✗ Branch 241 → 243 not taken.
✗ Branch 243 → 244 not taken.
✗ Branch 243 → 245 not taken.
✗ Branch 245 → 246 not taken.
✗ Branch 245 → 247 not taken.
✗ Branch 251 → 252 not taken.
✗ Branch 251 → 253 not taken.
✗ Branch 253 → 254 not taken.
✗ Branch 253 → 255 not taken.
✗ Branch 255 → 256 not taken.
✗ Branch 255 → 257 not taken.
✗ Branch 257 → 258 not taken.
✗ Branch 257 → 259 not taken.
✗ Branch 259 → 260 not taken.
✗ Branch 259 → 261 not taken.
✗ Branch 261 → 262 not taken.
✗ Branch 261 → 274 not taken.
✗ Branch 263 → 264 not taken.
✗ Branch 263 → 273 not taken.
✗ Branch 264 → 265 not taken.
✗ Branch 264 → 266 not taken.
✗ Branch 266 → 267 not taken.
✗ Branch 266 → 268 not taken.
✗ Branch 268 → 269 not taken.
✗ Branch 268 → 270 not taken.
✗ Branch 270 → 271 not taken.
✗ Branch 270 → 272 not taken.
✓ Branch 275 → 276 taken 80 times.
✓ Branch 275 → 277 taken 38 times.
✓ Branch 277 → 278 taken 80 times.
✓ Branch 277 → 279 taken 38 times.
✓ Branch 284 → 285 taken 80 times.
✓ Branch 284 → 286 taken 38 times.
✓ Branch 286 → 287 taken 80 times.
✓ Branch 286 → 299 taken 38 times.
✓ Branch 288 → 289 taken 240 times.
✓ Branch 288 → 298 taken 80 times.
✓ Branch 289 → 290 taken 240 times.
✗ Branch 289 → 291 not taken.
✓ Branch 291 → 292 taken 240 times.
✗ Branch 291 → 293 not taken.
✓ Branch 293 → 294 taken 240 times.
✗ Branch 293 → 295 not taken.
✓ Branch 295 → 296 taken 240 times.
✗ Branch 295 → 297 not taken.
✗ Branch 299 → 300 not taken.
✓ Branch 299 → 310 taken 118 times.
✗ Branch 301 → 302 not taken.
✗ Branch 301 → 309 not taken.
✗ Branch 302 → 303 not taken.
✗ Branch 302 → 304 not taken.
✗ Branch 304 → 305 not taken.
✗ Branch 304 → 306 not taken.
✗ Branch 306 → 307 not taken.
✗ Branch 306 → 308 not taken.
✓ Branch 310 → 311 taken 80 times.
✓ Branch 310 → 321 taken 38 times.
✓ Branch 312 → 313 taken 240 times.
✓ Branch 312 → 320 taken 80 times.
✓ Branch 313 → 314 taken 240 times.
✗ Branch 313 → 315 not taken.
✓ Branch 315 → 316 taken 240 times.
✗ Branch 315 → 317 not taken.
✓ Branch 317 → 318 taken 240 times.
✗ Branch 317 → 319 not taken.
✓ Branch 321 → 322 taken 80 times.
✓ Branch 321 → 328 taken 38 times.
✓ Branch 323 → 324 taken 80 times.
✓ Branch 323 → 327 taken 80 times.
✓ Branch 324 → 325 taken 80 times.
✗ Branch 324 → 326 not taken.
✓ Branch 328 → 329 taken 80 times.
✓ Branch 328 → 335 taken 38 times.
✓ Branch 330 → 331 taken 80 times.
✓ Branch 330 → 334 taken 80 times.
✓ Branch 331 → 332 taken 80 times.
✗ Branch 331 → 333 not taken.
✓ Branch 335 → 336 taken 80 times.
✓ Branch 335 → 337 taken 38 times.
1154 end subroutine insert_polar_extractor_1forms_enrich
743
744 371 subroutine finalize_polar_extractor_1forms(this)
745 implicit none
746
747 class(MFormPolarRegularity), intent(inout) :: this
748
749 integer :: i, j, k, row_3d, s1, polar_index, col, nr_polar_splines
750 integer :: i_min(3), i_max(3), j_min(3), j_max(3), k_min(3), k_max(3), col0, nr_nonzeros(3)
751 real(wp) :: coeff
752
753
2/2
✓ Branch 3 → 4 taken 1113 times.
✓ Branch 3 → 6 taken 371 times.
1484 do s1 = 1, 3
754 1484 call this%index_list(s1)%get_bounds(i_min(s1), i_max(s1), j_min(s1), j_max(s1), k_min(s1), k_max(s1))
755 end do
756
757 371 nr_nonzeros(1) = size(this%index_list(1)) * size(this%E1, 4)
758 371 nr_nonzeros(2) = size(this%index_list(2)) * size(this%E1, 4)
759
2/2
✓ Branch 7 → 8 taken 241 times.
✓ Branch 7 → 9 taken 130 times.
371 if (this%space%nr_components == 2) then
760 241 nr_nonzeros(3) = 0
761 nr_polar_splines = size(this%E1, 4)
762 else
763 130 nr_nonzeros(3) = size(this%index_list(3)) * size(this%E2, 4)
764 130 nr_polar_splines = size(this%E1, 4) + size(this%E2, 4)
765 end if
766
2/2
✓ Branch 11 → 12 taken 1113 times.
✓ Branch 11 → 14 taken 371 times.
1484 do s1 = 1, 3
767 ! TODO is using k_max(1) OK with MPI? nr cols needs to be same for all blocks
768 call this%extractor(s1)%init(nr_rows=size(this%index_list(s1)), nr_cols=nr_polar_splines * &
769 1484 & size(this%index_list(1), dir=3), nr_nonzero=nr_nonzeros(s1))
770 end do
771
772
2/2
✓ Branch 15 → 16 taken 1113 times.
✓ Branch 15 → 45 taken 371 times.
1484 do s1 = 1, 3
773
4/4
✓ Branch 16 → 17 taken 723 times.
✓ Branch 16 → 18 taken 390 times.
✓ Branch 17 → 18 taken 482 times.
✓ Branch 17 → 43 taken 241 times.
1113 if (this%space%nr_components == 2 .and. s1 == 3) cycle
774
775 742 select case (s1)
776 case (1, 2)
777 742 nr_polar_splines = size(this%E1, 4)
778 742 col0 = 0
779 case (3)
780 130 nr_polar_splines = size(this%E2, 4)
781
2/2
✓ Branch 18 → 19 taken 742 times.
✓ Branch 18 → 20 taken 130 times.
872 col0 = size(this%E1, 4) * size(this%index_list(1), dir=3)
782 end select
783
784
2/2
✓ Branch 22 → 23 taken 872 times.
✓ Branch 22 → 24 taken 2558 times.
3801 do k = k_min(s1), k_max(s1)
785
2/2
✓ Branch 25 → 26 taken 67816 times.
✓ Branch 25 → 41 taken 2558 times.
71487 do j = j_min(s1), j_max(s1)
786
2/2
✓ Branch 27 → 28 taken 160228 times.
✓ Branch 27 → 39 taken 67816 times.
230602 do i = i_min(s1), i_max(s1)
787 160228 row_3d = this%index_list(s1)%get_local_index(i, j, k)
788 polar_index = 0
789 ! NOTE: we skip gamma == 0 because it's gradient is zero
790
2/2
✓ Branch 29 → 30 taken 1699032 times.
✓ Branch 29 → 37 taken 160228 times.
1927076 do polar_index = 0, nr_polar_splines - 1
791 1699032 col = polar_index + (k - k_min(s1)) * nr_polar_splines
792
793 1196128 select case (s1)
794 case (1, 2)
795 1196128 coeff = this%E1(i, j, s1, polar_index + 1)
796 case (3)
797
2/2
✓ Branch 30 → 31 taken 1196128 times.
✓ Branch 30 → 32 taken 502904 times.
1699032 coeff = this%E2(i, j, s1, polar_index + 1)
798 end select
799
800
2/2
✓ Branch 33 → 34 taken 1032299 times.
✓ Branch 33 → 36 taken 666733 times.
1859260 if (abs(coeff) > COEFF_THR) call this%extractor(s1)%insert(row_3d, col0 + col, coeff)
801 end do
802 end do
803 end do
804 end do
805 end do
806
807
2/2
✓ Branch 46 → 47 taken 1113 times.
✓ Branch 46 → 49 taken 371 times.
1484 do s1 = 1, 3
808 1484 call this%extractor(s1)%finalize()
809 end do
810 371 end subroutine finalize_polar_extractor_1forms
811
812 !> @brief Initialize the 2-form extractor for the polar regularity constraint with curls of polar one-forms
813 !>
814 !> @param[inout] this The MFormPolarRegularity object to initialize
815 249 subroutine insert_polar_extractor_2forms(this)
816 use m_tensorprod, only: size
817 use m_tensorprod_indices, only: TensorProdIndices
818 use m_qr_factorisation, only: dense_column_space
819 implicit none
820
821 class(MFormPolarRegularity), intent(inout) :: this
822
823 type(TensorProdIndices) :: inds
824
825 integer :: nr_nonzeros(3), nr_polar_splines(3), nr_polar_splines_xy, nr_polar_splines_z, nr_bspline_y
826 integer :: i, j, j_plus_one, s1, polar_index, col0, col, row_xy, row_xy0, non_zeros_xy
827 integer :: i_min(3), i_max(3), j_min(3), j_max(3), k_min(3), k_max(3)
828 integer :: gamma, alpha, gamma_min, gamma_max
829 real(wp) :: coeff
830
831 249 nr_bspline_y = size(this%space%tp_spaces(1), 2)
832
833 gamma_min = 1
834 249 gamma_max = this%regularity
835
836 !> Initialize the indices on which this extractor acts
837 249 i_max(1) = min(this%regularity, this%space%tp_spaces(1)%rank_resp_bspline%i1)
838 249 i_max(2) = min(this%regularity - 1, this%space%tp_spaces(2)%rank_resp_bspline%i1)
839 ! NOTE the 3rd component gets at least 0, because of the pushforward factor 1/xp
840 249 i_max(3) = min(this%regularity, this%space%tp_spaces(3)%spaces(1)%degree, this%space%tp_spaces(3)%rank_resp_bspline%i1)
841 249 this%row0_2d(1) = 0
842
2/2
✓ Branch 3 → 4 taken 747 times.
✓ Branch 3 → 7 taken 249 times.
996 do s1 = 1, 3
843 747 i_min(s1) = max(0, this%space%tp_spaces(s1)%rank_resp_bspline%i0)
844 747 j_min(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%j0; j_max(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%j1
845 747 k_min(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%k0; k_max(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%k1
846
847 747 call inds%init(i_min(s1), i_max(s1), j_min(s1), j_max(s1), k_min(s1), k_max(s1), this%space%tp_spaces(s1)%spaces)
848 747 call this%index_list(s1)%init(inds, this%space%tp_spaces(s1)%rank_resp_bspline, this%space%tp_spaces(s1)%rank_l0)
849 996 this%row0_2d(s1 + 1) = this%row0_2d(s1) + size_2d(this%index_list(s1))
850 end do
851
852 ! Assume no enrichment: can be overridden later
853 249 call this%init_shmem_window(this%E2, this%shmem_E2, 0)
854
855
2/2
✓ Branch 9 → 10 taken 125 times.
✓ Branch 9 → 12 taken 124 times.
249 if (this%space%nr_components == 2) then
856 125 call this%init_shmem_window(this%E1, this%shmem_E1, 0)
857 125 return
858 end if
859
860 !> Initialize the radial and angular B-splines (data members of this object, using space0!)
861 call this%compute_1d_bsplines(this%space0%tp_spaces(1)%spaces(1), this%space0%tp_spaces(1)%spaces(2), 0, this%regularity, &
862 124 radial_truncate=this%regularity)
863
864 ! Curls of the Z-component of polar one-forms; resulting in the xy-components of the 2-form space
865 124 nr_polar_splines_xy = get_nr_polar_zeroforms(gamma_max) - get_nr_polar_zeroforms(0)
866
867 124 call this%init_shmem_window(this%E1, this%shmem_E1, nr_polar_splines_xy)
868
2/2
✓ Branch 14 → 15 taken 248 times.
✓ Branch 14 → 36 taken 124 times.
372 do s1 = 1, 2
869
2/2
✓ Branch 16 → 17 taken 7400 times.
✓ Branch 16 → 34 taken 248 times.
7772 do j = j_min(s1), j_max(s1)
870
2/2
✓ Branch 18 → 19 taken 16900 times.
✓ Branch 18 → 32 taken 7400 times.
24548 do i = i_min(s1), i_max(s1)
871 polar_index = 0
872
2/2
✓ Branch 19 → 20 taken 47176 times.
✓ Branch 19 → 30 taken 16900 times.
71476 do gamma = gamma_min, gamma_max
873
1/2
✗ Branch 20 → 21 not taken.
✓ Branch 20 → 22 taken 47176 times.
64076 do alpha = -gamma, gamma, 2
874
875 87976 select case (s1)
876 case (1)
877 87976 j_plus_one = mod(j + 1, nr_bspline_y)
878 coeff = (this%radial_bsplines(gamma)%data(i) * &
879 87976 (this%angular_bsplines(alpha)%data(j_plus_one) - this%angular_bsplines(alpha)%data(j)))
880 case (2)
881 coeff = -(this%radial_bsplines(gamma)%data(i + 1) - this%radial_bsplines(gamma)%data(i)) * &
882
2/2
✓ Branch 22 → 23 taken 87976 times.
✓ Branch 22 → 24 taken 67932 times.
155908 this%angular_bsplines(alpha)%data(j)
883 end select
884
2/2
✓ Branch 25 → 26 taken 71593 times.
✓ Branch 25 → 27 taken 84315 times.
155908 if (abs(coeff) > COEFF_THR) this%E1(i, j, s1, polar_index + 1) = coeff
885
2/2
✓ Branch 27 → 22 taken 108732 times.
✓ Branch 27 → 28 taken 47176 times.
155908 polar_index = polar_index + 1
886 end do
887 end do
888 end do
889 end do
890 end do
891 end subroutine insert_polar_extractor_2forms
892
893 !> @brief Enrich 2-form basis for the polar regularity constraint using weak gradients of polar 3-forms
894 !>
895 !> @param[inout] this The MFormPolarEnriched object to initialize
896 106 subroutine insert_polar_extractor_2forms_enrich(this)
897 implicit none
898
899 class(MFormPolarRegularity), intent(inout) :: this
900
901 integer :: i, j, polar_index, alpha, gamma, gamma_min, gamma_max, nr_polar_splines_z
902 integer :: i_min(3), i_max(3), j_min(3), j_max(3), s1
903
904 ! TODO include weak gradients of 3-forms for enrichment
905
906 gamma_min = 0
907
908 ! The pushforward yields extra factor xp, and the 3rd component of the 2-form space has degree reduced by 1
909 106 gamma_max = this%regularity - 2
910
2/2
✓ Branch 2 → 3 taken 34 times.
✓ Branch 2 → 4 taken 72 times.
106 if (gamma_max < gamma_min) return
911
912 72 nr_polar_splines_z = get_nr_polar_zeroforms(gamma_max)
913
914
2/2
✓ Branch 5 → 6 taken 216 times.
✓ Branch 5 → 8 taken 72 times.
288 do s1 = 1, 3
915 288 call this%index_list(s1)%get_bounds(i_min(s1), i_max(s1), j_min(s1), j_max(s1))
916 end do
917
918 !> Initialize the radial and angular B-splines (data members of this object)
919 call this%compute_1d_bsplines(this%space%tp_spaces(3)%spaces(1), this%space%tp_spaces(3)%spaces(2), &
920 72 gamma_min + 1, gamma_max + 1, radial_truncate=gamma_max + 1)
921
922 72 call this%init_shmem_window(this%E2, this%shmem_E2, nr_polar_splines_z)
923 s1 = 3
924
2/2
✓ Branch 12 → 13 taken 2004 times.
✓ Branch 12 → 25 taken 72 times.
2076 do j = j_min(s1), j_max(s1)
925
2/2
✓ Branch 14 → 15 taken 7184 times.
✓ Branch 14 → 23 taken 2004 times.
9260 do i = i_min(s1), i_max(s1)
926 polar_index = 0
927
2/2
✓ Branch 15 → 16 taken 16980 times.
✓ Branch 15 → 21 taken 7184 times.
26168 do gamma = gamma_min, gamma_max
928
1/2
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 16980 times.
24164 do alpha = -gamma, gamma, 2
929 this%E2(i, j, s1, polar_index + 1) = &
930
2/2
✓ Branch 18 → 18 taken 16538 times.
✓ Branch 18 → 19 taken 16980 times.
33518 this%radial_bsplines(gamma + 1)%data(i) * this%angular_bsplines(alpha)%data(j)
931 16980 polar_index = polar_index + 1
932 end do
933 end do
934 end do
935 end do
936 end subroutine insert_polar_extractor_2forms_enrich
937
938 249 subroutine finalize_polar_extractor_2forms(this)
939 implicit none
940
941 class(MFormPolarRegularity), intent(inout) :: this
942
943 integer :: s1
944 integer :: i, j, k, row_3d, polar_index, col, nr_polar_splines
945 integer :: i_min(3), i_max(3), j_min(3), j_max(3), k_min(3), k_max(3), col0, nr_nonzeros(3)
946 real(wp) :: coeff
947
948
2/2
✓ Branch 3 → 4 taken 747 times.
✓ Branch 3 → 6 taken 249 times.
996 do s1 = 1, 3
949 996 call this%index_list(s1)%get_bounds(i_min(s1), i_max(s1), j_min(s1), j_max(s1), k_min(s1), k_max(s1))
950 end do
951
952 249 nr_nonzeros(1) = size(this%index_list(1)) * size(this%E1, 4)
953 249 nr_nonzeros(2) = size(this%index_list(2)) * size(this%E1, 4)
954 249 nr_nonzeros(3) = size(this%index_list(3)) * size(this%E2, 4)
955 249 nr_polar_splines = size(this%E1, 4) + size(this%E2, 4)
956
957
2/2
✓ Branch 8 → 9 taken 747 times.
✓ Branch 8 → 11 taken 249 times.
996 do s1 = 1, 3
958 call this%extractor(s1)%init(nr_rows=size(this%index_list(s1)), nr_cols=nr_polar_splines * &
959 996 & size(this%index_list(3), dir=3), nr_nonzero=nr_nonzeros(s1))
960 end do
961
962
2/2
✓ Branch 12 → 13 taken 747 times.
✓ Branch 12 → 41 taken 249 times.
996 do s1 = 1, 3
963
2/2
✓ Branch 13 → 14 taken 498 times.
✓ Branch 13 → 16 taken 249 times.
747 if (s1 < 3) then
964
2/2
✓ Branch 14 → 15 taken 248 times.
✓ Branch 14 → 39 taken 250 times.
498 if (this%space%nr_components == 2) cycle
965 248 nr_polar_splines = size(this%E1, 4)
966 col0 = 0
967 else
968 249 nr_polar_splines = size(this%E2, 4)
969 249 col0 = size(this%E1, 4) * size(this%index_list(1), dir=3)
970 end if
971
972
2/2
✓ Branch 18 → 19 taken 497 times.
✓ Branch 18 → 20 taken 1985 times.
2731 do k = k_min(s1), k_max(s1)
973
2/2
✓ Branch 21 → 22 taken 54230 times.
✓ Branch 21 → 37 taken 1985 times.
56962 do j = j_min(s1), j_max(s1)
974
2/2
✓ Branch 23 → 24 taken 133022 times.
✓ Branch 23 → 35 taken 54230 times.
189237 do i = i_min(s1), i_max(s1)
975 133022 row_3d = this%index_list(s1)%get_local_index(i, j, k)
976 polar_index = 0
977
978
2/2
✓ Branch 25 → 26 taken 859266 times.
✓ Branch 25 → 33 taken 133022 times.
1046518 do polar_index = 0, nr_polar_splines - 1
979 859266 col = polar_index + (k - k_min(s1)) * nr_polar_splines
980
2/2
✓ Branch 26 → 27 taken 756584 times.
✓ Branch 26 → 28 taken 102682 times.
859266 if (s1 < 3) then
981 756584 coeff = this%E1(i, j, s1, polar_index + 1)
982 else
983 102682 coeff = this%E2(i, j, s1, polar_index + 1)
984 end if
985
986
2/2
✓ Branch 29 → 30 taken 428183 times.
✓ Branch 29 → 32 taken 431083 times.
992288 if (abs(coeff) > COEFF_THR) call this%extractor(s1)%insert(row_3d, col0 + col, coeff)
987 end do
988 end do
989 end do
990 end do
991 end do
992
993
2/2
✓ Branch 42 → 43 taken 747 times.
✓ Branch 42 → 45 taken 249 times.
996 do s1 = 1, 3
994 996 call this%extractor(s1)%finalize()
995 end do
996
997 249 end subroutine finalize_polar_extractor_2forms
998
999 !> @brief Initialize the 3-form extractor for the polar regularity constraint
1000 !>
1001 !> @param[inout] this The MFormPolarRegularity object to initialize
1002 124 subroutine insert_polar_extractor_3forms(this)
1003 use m_tensorprod, only: size
1004 use m_tensorprod_indices, only: TensorProdIndices
1005 implicit none
1006
1007 class(MFormPolarRegularity), intent(inout) :: this
1008
1009 type(TensorProdIndices) :: inds
1010
1011 integer :: nr_bspline_y, i_min, i_max, j_min, j_max, k_min, k_max
1012
1013 124 nr_bspline_y = size(this%space%tp_spaces(1), 2)
1014
1015 !> Initialize the indices on which this extractor acts
1016 124 i_min = max(0, this%space%tp_spaces(1)%rank_resp_bspline%i0)
1017 ! NOTE i_max is at least 0, because of the pushforward factor 1/xp
1018 124 i_max = min(this%regularity, this%space%tp_spaces(1)%spaces(1)%degree, this%space%tp_spaces(1)%rank_resp_bspline%i1)
1019 124 j_min = this%space%tp_spaces(1)%rank_resp_bspline%j0; j_max = this%space%tp_spaces(1)%rank_resp_bspline%j1
1020 124 k_min = this%space%tp_spaces(1)%rank_resp_bspline%k0; k_max = this%space%tp_spaces(1)%rank_resp_bspline%k1
1021
1022 124 call inds%init(i_min, i_max, j_min, j_max, k_min, k_max, this%space%tp_spaces(1)%spaces)
1023 124 call this%index_list(1)%init(inds, this%space%tp_spaces(1)%rank_resp_bspline, this%space%tp_spaces(1)%rank_l0)
1024
1025 124 call this%init_shmem_window(this%E1, this%shmem_E1, 0)
1026 124 end subroutine insert_polar_extractor_3forms
1027
1028 !> @brief Enrich 3-form basis for the polar regularity constraint using regular 3-forms
1029 !>
1030 !> @param[inout] this The MFormPolarEnriched object to initialize
1031 55 subroutine insert_polar_extractor_3forms_enrich(this)
1032 implicit none
1033
1034 class(MFormPolarRegularity), intent(inout) :: this
1035
1036 integer :: i, j, polar_index, alpha, gamma, gamma_min, gamma_max, nr_polar_splines_z
1037 integer :: i_min, i_max, j_min, j_max
1038
1039 gamma_min = 0
1040
1041 ! The pushforward yields extra factor xp, and the 3rd component of the 2-form space has degree reduced by 1
1042 55 gamma_max = this%regularity - 2
1043
2/2
✓ Branch 2 → 3 taken 18 times.
✓ Branch 2 → 4 taken 37 times.
55 if (gamma_max < gamma_min) return
1044
1045 37 nr_polar_splines_z = get_nr_polar_zeroforms(gamma_max)
1046
1047 37 call this%index_list(1)%get_bounds(i_min, i_max, j_min, j_max)
1048
1049 !> Initialize the radial and angular B-splines (data members of this object)
1050 call this%compute_1d_bsplines(this%space%tp_spaces(1)%spaces(1), this%space%tp_spaces(1)%spaces(2), &
1051 37 gamma_min + 1, gamma_max + 1, radial_truncate=gamma_max + 1)
1052
1053 37 call this%init_shmem_window(this%E1, this%shmem_E1, nr_polar_splines_z)
1054
1055
2/2
✓ Branch 8 → 9 taken 1008 times.
✓ Branch 8 → 21 taken 37 times.
1045 do j = j_min, j_max
1056
2/2
✓ Branch 10 → 11 taken 3556 times.
✓ Branch 10 → 19 taken 1008 times.
4601 do i = i_min, i_max
1057 polar_index = 0
1058
2/2
✓ Branch 11 → 12 taken 8184 times.
✓ Branch 11 → 17 taken 3556 times.
12748 do gamma = gamma_min, gamma_max
1059
1/2
✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 8184 times.
11740 do alpha = -gamma, gamma, 2
1060 this%E1(i, j, 1, polar_index + 1) = &
1061
2/2
✓ Branch 14 → 14 taken 7684 times.
✓ Branch 14 → 15 taken 8184 times.
15868 this%radial_bsplines(gamma + 1)%data(i) * this%angular_bsplines(alpha)%data(j)
1062 8184 polar_index = polar_index + 1
1063 end do
1064 end do
1065 end do
1066 end do
1067 end subroutine insert_polar_extractor_3forms_enrich
1068
1069 124 subroutine finalize_polar_extractor_3forms(this)
1070 implicit none
1071
1072 class(MFormPolarRegularity), intent(inout) :: this
1073
1074 integer :: i, j, k, row_2d, row_3d, polar_index, col, nr_polar_splines
1075 integer :: i_min, i_max, j_min, j_max, k_min, k_max, nr_nonzeros
1076 real(wp) :: coeff
1077
1078 124 call this%index_list(1)%get_bounds(i_min, i_max, j_min, j_max, k_min, k_max)
1079
1080 124 nr_nonzeros = size(this%index_list(1)) * size(this%E1, 4)
1081 nr_polar_splines = size(this%E1, 4)
1082
1083 call this%extractor(1)%init(nr_rows=size(this%index_list(1)), nr_cols=nr_polar_splines * &
1084 124 & size(this%index_list(1), dir=3), nr_nonzero=nr_nonzeros)
1085
1086
2/2
✓ Branch 5 → 6 taken 620 times.
✓ Branch 5 → 22 taken 124 times.
744 do k = k_min, k_max
1087
2/2
✓ Branch 7 → 8 taken 16800 times.
✓ Branch 7 → 20 taken 620 times.
17544 do j = j_min, j_max
1088
2/2
✓ Branch 9 → 10 taken 42592 times.
✓ Branch 9 → 18 taken 16800 times.
60012 do i = i_min, i_max
1089 42592 row_3d = this%index_list(1)%get_local_index(i, j, k)
1090 row_2d = this%index_list(1)%get_local_index(i, j)
1091 polar_index = 0
1092
1093
2/2
✓ Branch 12 → 13 taken 85032 times.
✓ Branch 12 → 16 taken 42592 times.
144424 do polar_index = 0, nr_polar_splines - 1
1094 85032 col = polar_index + (k - k_min) * nr_polar_splines
1095 85032 coeff = this%E1(i, j, 1, polar_index + 1)
1096
1097
2/2
✓ Branch 13 → 11 taken 47244 times.
✓ Branch 13 → 14 taken 37788 times.
127624 if (abs(coeff) > COEFF_THR) call this%extractor(1)%insert(row_3d, col, coeff)
1098 end do
1099 end do
1100 end do
1101 end do
1102
1103 124 call this%extractor(1)%finalize()
1104
1105 124 end subroutine finalize_polar_extractor_3forms
1106
1107 !> @brief Compute the 1D B-spline coefficients for the polar regularity constraint
1108 !>
1109 !> @param[inout] this The MFormPolarRegularity object for which to compute the B-splines
1110 !> @param[in] bspline_x The B-spline space in the x-direction
1111 !> @param[in] bspline_y The B-spline space in the y-direction
1112 !> @param[in] gamma_min The minimum radial degree
1113 !> @param[in] gamma_max The maximum radial degree
1114 !> @param[in] radial_truncate _(optional)_ The index at which to truncate the radial functions, 0:radial_truncate is kept
1115 1217 subroutine compute_1d_bsplines(this, bspline_x, bspline_y, gamma_min, gamma_max, radial_truncate)
1116 use m_bspline, only: BSplineSpace, size
1117 use m_bspline_linalg, only: l2_projection, biorthogonalize
1118 use m_common, only: user_function_1d_interface
1119
1120 implicit none
1121
1122 class(MFormPolarRegularity), intent(inout) :: this
1123 type(BSplineSpace), intent(in) :: bspline_x, bspline_y
1124 integer, intent(in) :: gamma_min, gamma_max
1125 integer, optional, intent(in) :: radial_truncate
1126
1127 integer :: alpha, gamma, nr_bsplines_y
1128
1129
2/4
✓ Branch 2 → 3 taken 1217 times.
✗ Branch 2 → 4 not taken.
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 1217 times.
1217 if (gamma_min < 0 .or. gamma_max > bspline_x%degree) then
1130 error stop "MFormPolarRegularity::compute_1d_bsplines: regularity must be in [0, degree_x] where degree_x is "// &
1131 & " the degree of the B-spline basis in the x-direction"
1132 end if
1133
1134
1/2
✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 1217 times.
1217 if (2 * gamma_max + 1 >= size(bspline_y)) then
1135 error stop "MFormPolarRegularity::compute_1d_bsplines: 2*regularity + 1 must be less than the number of basis"// &
1136 & " functions in the y-direction"
1137 end if
1138
1139
6/8
✓ Branch 7 → 8 taken 117 times.
✓ Branch 7 → 16 taken 1100 times.
✓ Branch 9 → 10 taken 467 times.
✓ Branch 9 → 13 taken 117 times.
✓ Branch 10 → 11 taken 467 times.
✗ Branch 10 → 12 not taken.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 15 taken 117 times.
1684 if (allocated(this%radial_bsplines)) deallocate (this%radial_bsplines)
1140
6/8
✓ Branch 16 → 17 taken 117 times.
✓ Branch 16 → 25 taken 1100 times.
✓ Branch 18 → 19 taken 817 times.
✓ Branch 18 → 22 taken 117 times.
✓ Branch 19 → 20 taken 817 times.
✗ Branch 19 → 21 not taken.
✗ Branch 22 → 23 not taken.
✓ Branch 22 → 24 taken 117 times.
2034 if (allocated(this%angular_bsplines)) deallocate (this%angular_bsplines)
1141
1142
9/16
✗ Branch 25 → 26 not taken.
✓ Branch 25 → 27 taken 1217 times.
✓ Branch 27 → 26 taken 1217 times.
✗ Branch 27 → 28 not taken.
✓ Branch 28 → 29 taken 1217 times.
✗ Branch 28 → 30 not taken.
✓ Branch 30 → 31 taken 1217 times.
✗ Branch 30 → 32 not taken.
✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 1217 times.
✗ Branch 34 → 35 not taken.
✓ Branch 34 → 36 taken 1217 times.
✗ Branch 36 → 37 not taken.
✓ Branch 36 → 38 taken 1217 times.
✓ Branch 39 → 40 taken 3089 times.
✓ Branch 39 → 41 taken 1217 times.
6740 allocate (this%radial_bsplines(gamma_min:gamma_max))
1143
9/16
✗ Branch 41 → 42 not taken.
✓ Branch 41 → 43 taken 1217 times.
✓ Branch 43 → 42 taken 1217 times.
✗ Branch 43 → 44 not taken.
✓ Branch 44 → 45 taken 1217 times.
✗ Branch 44 → 46 not taken.
✓ Branch 46 → 47 taken 1217 times.
✗ Branch 46 → 48 not taken.
✗ Branch 48 → 49 not taken.
✓ Branch 48 → 50 taken 1217 times.
✗ Branch 50 → 51 not taken.
✓ Branch 50 → 52 taken 1217 times.
✗ Branch 52 → 53 not taken.
✓ Branch 52 → 54 taken 1217 times.
✓ Branch 55 → 56 taken 5339 times.
✓ Branch 55 → 57 taken 1217 times.
8990 allocate (this%angular_bsplines(-gamma_max:gamma_max))
1144
1145 ! Compute L2 approximation of the radial and angular functions which make up the Zernike basis
1146
2/2
✓ Branch 58 → 59 taken 3089 times.
✓ Branch 58 → 84 taken 1217 times.
4306 do gamma = gamma_min, gamma_max
1147
2/2
✓ Branch 59 → 60 taken 1028 times.
✓ Branch 59 → 64 taken 2061 times.
3089 if (gamma == 0) then
1148 1028 call this%radial_bsplines(gamma)%init(bspline_x)
1149
2/2
✓ Branch 62 → 63 taken 21420 times.
✓ Branch 62 → 68 taken 1028 times.
22448 this%radial_bsplines(gamma)%data(0:) = 1._wp
1150 else
1151 2061 call l2_projection(this%radial_bsplines(gamma), bspline_x, radial_fun)
1152 ! Analytically the first gamma coefficients are zero
1153
2/2
✓ Branch 66 → 67 taken 4478 times.
✓ Branch 66 → 68 taken 2061 times.
6539 this%radial_bsplines(gamma)%data(0:gamma - 1) = 0._wp
1154 end if
1155
3/4
✓ Branch 68 → 69 taken 3089 times.
✗ Branch 68 → 72 not taken.
✓ Branch 70 → 71 taken 51902 times.
✓ Branch 70 → 72 taken 3089 times.
54991 if (present(radial_truncate)) this%radial_bsplines(gamma)%data(radial_truncate + 1:) = 0._wp
1156
1157 ! Optional: we normalize the coefficients for numerical stability (merely a change of basis)
1158 this%radial_bsplines(gamma)%data(0:) = this%radial_bsplines(gamma)%data(0:) / &
1159
6/10
✓ Branch 73 → 74 taken 3089 times.
✗ Branch 73 → 76 not taken.
✗ Branch 74 → 75 not taken.
✓ Branch 74 → 78 taken 3089 times.
✗ Branch 76 → 77 not taken.
✗ Branch 76 → 78 not taken.
✓ Branch 79 → 80 taken 63566 times.
✓ Branch 79 → 81 taken 3089 times.
✓ Branch 81 → 82 taken 63566 times.
✓ Branch 81 → 83 taken 3089 times.
134527 maxval(abs(this%radial_bsplines(gamma)%data(0:)))
1160 end do
1161
1162 1217 nr_bsplines_y = size(bspline_y)
1163
2/2
✓ Branch 86 → 87 taken 5339 times.
✓ Branch 86 → 100 taken 1217 times.
6556 do alpha = -gamma_max, gamma_max
1164
2/2
✓ Branch 87 → 88 taken 1217 times.
✓ Branch 87 → 92 taken 4122 times.
6556 if (alpha == 0) then
1165 1217 call this%angular_bsplines(alpha)%init(bspline_y)
1166
2/2
✓ Branch 90 → 91 taken 39825 times.
✓ Branch 90 → 99 taken 1217 times.
41042 this%angular_bsplines(alpha)%data(:) = 1._wp
1167 else
1168 4122 call l2_projection(this%angular_bsplines(alpha), bspline_y, angular_fun)
1169 ! Force zero mean
1170 this%angular_bsplines(alpha)%data(:) = this%angular_bsplines(alpha)%data(:) - &
1171
4/4
✓ Branch 94 → 95 taken 120096 times.
✓ Branch 94 → 96 taken 4122 times.
✓ Branch 97 → 98 taken 134950 times.
✓ Branch 97 → 99 taken 4122 times.
259168 sum(this%angular_bsplines(alpha)%data(0:nr_bsplines_y - 1)) / nr_bsplines_y
1172 end if
1173 end do
1174
1175 contains
1176 177512 pure real(wp) function radial_fun(xp) result(ans)
1177 implicit none
1178 real(wp), intent(in) :: xp
1179
1180 177512 ans = xp**gamma
1181 177512 end function radial_fun
1182
1183 565792 pure real(wp) function angular_fun(yp) result(ans)
1184 use m_common, only: pi
1185
1186 implicit none
1187 real(wp), intent(in) :: yp
1188
1189
1/2
✓ Branch 2 → 3 taken 565792 times.
✗ Branch 2 → 6 not taken.
565792 if (alpha == 0) then
1190 ans = 1._wp
1191
2/2
✓ Branch 3 → 4 taken 282896 times.
✓ Branch 3 → 5 taken 282896 times.
565792 else if (alpha < 0) then
1192 282896 ans = sin(-2 * pi * alpha * yp)
1193 else
1194 282896 ans = cos(2 * pi * alpha * yp)
1195 end if
1196
1197 565792 end function angular_fun
1198 end subroutine compute_1d_bsplines
1199
1200 !> @brief Destroy the MFormPolarRegularity object
1201 !>
1202 !> @param[inout] this The MFormPolarRegularity object to destroy
1203 1766 subroutine destroy_mform_polarregularity(this)
1204 implicit none
1205 class(MFormPolarRegularity), intent(inout) :: this
1206 integer :: s1, regularity
1207
1208
2/2
✓ Branch 2 → 3 taken 523 times.
✓ Branch 2 → 21 taken 1243 times.
1766 if (allocated(this%radial_bsplines)) then
1209
4/6
✓ Branch 3 → 4 taken 523 times.
✗ Branch 3 → 5 not taken.
✓ Branch 5 → 6 taken 523 times.
✗ Branch 5 → 7 not taken.
✓ Branch 8 → 9 taken 1407 times.
✓ Branch 8 → 11 taken 523 times.
2976 do s1 = lbound(this%radial_bsplines, 1), ubound(this%radial_bsplines, 1)
1210 1930 call this%radial_bsplines(s1)%destroy()
1211 end do
1212
5/8
✓ Branch 12 → 13 taken 523 times.
✗ Branch 12 → 18 not taken.
✓ Branch 14 → 15 taken 1407 times.
✓ Branch 14 → 18 taken 523 times.
✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 1407 times.
✗ Branch 18 → 19 not taken.
✓ Branch 18 → 20 taken 523 times.
1930 deallocate (this%radial_bsplines)
1213 end if
1214
1215
2/2
✓ Branch 21 → 22 taken 523 times.
✓ Branch 21 → 40 taken 1243 times.
1766 if (allocated(this%angular_bsplines)) then
1216
4/6
✓ Branch 22 → 23 taken 523 times.
✗ Branch 22 → 24 not taken.
✓ Branch 24 → 25 taken 523 times.
✗ Branch 24 → 26 not taken.
✓ Branch 27 → 28 taken 2399 times.
✓ Branch 27 → 30 taken 523 times.
3968 do s1 = lbound(this%angular_bsplines, 1), ubound(this%angular_bsplines, 1)
1217 2922 call this%angular_bsplines(s1)%destroy()
1218 end do
1219
5/8
✓ Branch 31 → 32 taken 523 times.
✗ Branch 31 → 37 not taken.
✓ Branch 33 → 34 taken 2399 times.
✓ Branch 33 → 37 taken 523 times.
✗ Branch 34 → 35 not taken.
✓ Branch 34 → 36 taken 2399 times.
✗ Branch 37 → 38 not taken.
✓ Branch 37 → 39 taken 523 times.
2922 deallocate (this%angular_bsplines)
1220 end if
1221
1222 1766 call this%shmem_E1%destroy()
1223 1766 call this%shmem_E2%destroy()
1224
2/2
✓ Branch 42 → 43 taken 1038 times.
✓ Branch 42 → 44 taken 728 times.
1766 if (associated(this%E1)) nullify (this%E1)
1225
2/2
✓ Branch 44 → 45 taken 1099 times.
✓ Branch 44 → 46 taken 667 times.
1766 if (associated(this%E2)) nullify (this%E2)
1226
1227 1766 this%my_rank_window = -1
1228 1766 this%comm_window = -1
1229 1766 end subroutine destroy_mform_polarregularity
1230
1231 pure integer function get_nr_polar_zeroforms(regularity) result(nr)
1232 implicit none
1233 integer, intent(in) :: regularity
1234
1235
3/4
✓ Branch 2 → 3 taken 241 times.
✓ Branch 2 → 4 taken 130 times.
✗ Branch 115 → 116 not taken.
✓ Branch 115 → 117 taken 80 times.
1217 nr = (regularity + 1) * (regularity + 2) / 2
1236 end function get_nr_polar_zeroforms
1237
1238
74/122
__m_mform_constraint_polar_MOD___copy_m_mform_constraint_polar_Mformpolarregularity:
✓ Branch 2 → 3 taken 447 times.
✗ Branch 2 → 70 not taken.
✓ Branch 3 → 4 taken 132 times.
✓ Branch 3 → 5 taken 315 times.
✓ Branch 6 → 7 taken 132 times.
✓ Branch 6 → 22 taken 315 times.
✓ Branch 8 → 9 taken 132 times.
✓ Branch 8 → 23 taken 132 times.
✓ Branch 9 → 10 taken 132 times.
✗ Branch 9 → 11 not taken.
✓ Branch 12 → 13 taken 132 times.
✗ Branch 12 → 14 not taken.
✓ Branch 15 → 16 taken 132 times.
✗ Branch 15 → 17 not taken.
✓ Branch 18 → 19 taken 132 times.
✗ Branch 18 → 20 not taken.
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 36 taken 447 times.
✗ Branch 25 → 26 not taken.
✗ Branch 25 → 37 not taken.
✗ Branch 26 → 27 not taken.
✗ Branch 26 → 28 not taken.
✗ Branch 29 → 30 not taken.
✗ Branch 29 → 31 not taken.
✗ Branch 32 → 33 not taken.
✗ Branch 32 → 34 not taken.
✓ Branch 37 → 38 taken 132 times.
✓ Branch 37 → 50 taken 315 times.
✓ Branch 39 → 40 taken 132 times.
✓ Branch 39 → 51 taken 132 times.
✓ Branch 40 → 41 taken 132 times.
✗ Branch 40 → 42 not taken.
✓ Branch 43 → 44 taken 132 times.
✗ Branch 43 → 45 not taken.
✓ Branch 46 → 47 taken 132 times.
✗ Branch 46 → 48 not taken.
✓ Branch 51 → 52 taken 132 times.
✓ Branch 51 → 58 taken 315 times.
✓ Branch 53 → 54 taken 492 times.
✓ Branch 53 → 59 taken 132 times.
✓ Branch 54 → 55 taken 492 times.
✗ Branch 54 → 56 not taken.
✓ Branch 59 → 60 taken 132 times.
✓ Branch 59 → 66 taken 315 times.
✓ Branch 61 → 62 taken 852 times.
✓ Branch 61 → 67 taken 132 times.
✓ Branch 62 → 63 taken 852 times.
✗ Branch 62 → 64 not taken.
✓ Branch 67 → 68 taken 132 times.
✓ Branch 67 → 69 taken 315 times.
__m_mform_constraint_polar_MOD___final_m_mform_constraint_polar_Mformpolarregularity:
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 7 taken 447 times.
✗ Branch 4 → 5 not taken.
✗ Branch 4 → 6 not taken.
✓ Branch 8 → 9 taken 447 times.
✓ Branch 8 → 76 taken 447 times.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 447 times.
✓ Branch 12 → 13 taken 447 times.
✗ Branch 12 → 15 not taken.
✓ Branch 13 → 14 taken 447 times.
✗ Branch 13 → 16 not taken.
✗ Branch 15 → 16 not taken.
✗ Branch 15 → 55 not taken.
✓ Branch 16 → 17 taken 2 times.
✓ Branch 16 → 28 taken 445 times.
✓ Branch 18 → 19 taken 2 times.
✓ Branch 18 → 28 taken 2 times.
✓ Branch 19 → 20 taken 2 times.
✗ Branch 19 → 21 not taken.
✓ Branch 21 → 22 taken 2 times.
✗ Branch 21 → 23 not taken.
✓ Branch 23 → 24 taken 2 times.
✗ Branch 23 → 25 not taken.
✓ Branch 25 → 26 taken 2 times.
✗ Branch 25 → 27 not taken.
✓ Branch 28 → 29 taken 2 times.
✓ Branch 28 → 30 taken 445 times.
✗ Branch 30 → 31 not taken.
✓ Branch 30 → 40 taken 447 times.
✗ Branch 32 → 33 not taken.
✗ Branch 32 → 40 not taken.
✗ Branch 33 → 34 not taken.
✗ Branch 33 → 35 not taken.
✗ Branch 35 → 36 not taken.
✗ Branch 35 → 37 not taken.
✗ Branch 37 → 38 not taken.
✗ Branch 37 → 39 not taken.
✗ Branch 40 → 41 not taken.
✓ Branch 40 → 42 taken 447 times.
✓ Branch 42 → 43 taken 2 times.
✓ Branch 42 → 52 taken 445 times.
✓ Branch 44 → 45 taken 2 times.
✓ Branch 44 → 52 taken 2 times.
✓ Branch 45 → 46 taken 2 times.
✗ Branch 45 → 47 not taken.
✓ Branch 47 → 48 taken 2 times.
✗ Branch 47 → 49 not taken.
✓ Branch 49 → 50 taken 2 times.
✗ Branch 49 → 51 not taken.
✓ Branch 52 → 53 taken 2 times.
✓ Branch 52 → 54 taken 445 times.
✓ Branch 55 → 56 taken 447 times.
✗ Branch 55 → 73 not taken.
✓ Branch 56 → 57 taken 1 time.
✓ Branch 56 → 62 taken 446 times.
✓ Branch 58 → 59 taken 9 times.
✓ Branch 58 → 62 taken 1 time.
✓ Branch 59 → 60 taken 9 times.
✗ Branch 59 → 61 not taken.
✓ Branch 62 → 63 taken 1 time.
✓ Branch 62 → 64 taken 446 times.
✓ Branch 64 → 65 taken 1 time.
✓ Branch 64 → 70 taken 446 times.
✓ Branch 66 → 67 taken 17 times.
✓ Branch 66 → 70 taken 1 time.
✓ Branch 67 → 68 taken 17 times.
✗ Branch 67 → 69 not taken.
✓ Branch 70 → 71 taken 1 time.
✓ Branch 70 → 72 taken 446 times.
✓ Branch 73 → 74 taken 447 times.
✗ Branch 73 → 75 not taken.
3426 end module m_mform_constraint_polar
1239