1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
use crate::errors::SGError;
use crate::utilities::float::Float;
use crate::{basis::{base::Basis, linear::LinearBasis}, dynamic::{iterators::dynamic_grid_iterator::GridIteratorT}};
#[inline]
#[allow(clippy::too_many_arguments)]
pub(crate) fn eval_boundary<T: Float +std::ops::AddAssign, Iterator: GridIteratorT>(x: &[f64],
dim: usize, value: T, iterator: &mut Iterator, alpha: &[T], result: &mut [T], ndim: usize, num_outputs: usize) -> Result<(), SGError>
{
// LinearBasis is a zero-sized type, create it inline instead of passing as parameter
let basis = LinearBasis;
let mut level = 0;
loop
{
let node_index = iterator.index().ok_or(SGError::InvalidIndex)?;
let work_index = iterator.point_index(dim);
if level > 0
{
let new_value = T::from(basis.eval(level, work_index, x[dim]));
if dim == ndim - 1
{
#[allow(clippy::needless_range_loop)]
for i in 0..num_outputs
{
result[i] += alpha[node_index*num_outputs+i] * value * new_value;
}
}
else
{
eval_boundary(x, dim + 1, value * new_value, iterator, alpha, result, ndim, num_outputs)?;
}
}
else
{
// handle boundaries if we are on level 0
// level 0, index 0
// reset_to_left_level_zero now checks if the node exists - after grid coarsening some boundary nodes are removed.
if iterator.reset_to_left_level_zero(dim)
{
let seq_l = iterator.index().ok_or(SGError::InvalidIndex)?;
let new_value_l = T::from(basis.eval(0, 0, x[dim]));
if dim == ndim - 1
{
#[allow(clippy::needless_range_loop)]
for i in 0..num_outputs
{
result[i] += alpha[seq_l*num_outputs+i] * value * new_value_l;
}
}
else
{
eval_boundary(x, dim + 1, value * new_value_l, iterator, alpha, result, ndim, num_outputs)?;
}
}
// reset_to_right_level_zero now checks if the node exists - after grid coarsening some boundary nodes are removed.
if iterator.reset_to_right_level_zero(dim)
{
let seq_r = iterator.index().ok_or(SGError::InvalidIndex)?;
let new_value_r = T::from(basis.eval(0, 1, x[dim]));
if dim == ndim - 1
{
#[allow(clippy::needless_range_loop)]
for i in 0..num_outputs
{
result[i] += alpha[seq_r*num_outputs+i] * value * new_value_r;
}
}
else
{
eval_boundary(x, dim + 1, value * new_value_r, iterator, alpha, result, ndim, num_outputs)?;
}
}
}
// Can't descend any further.
if iterator.is_leaf()
{
break;
}
// this decides in which direction we should descend by evaluating
// the corresponding bit
// the bits are coded from left to right starting with level 1
// being in position max_level
let x = x[dim];
if level > 0
{
let h = 1 << level;
let x_coord = work_index as f64 / h as f64;
if (x - x_coord).abs() < 1e-15
{
break;
}
if x > x_coord && !iterator.right_child(dim)
{
break;
}
if x <= x_coord && !iterator.left_child(dim)
{
break;
}
}
else
{
if x.abs() < 1e-15 || (x-1.0).abs() < 1e-15
{
break;
}
if !iterator.reset_to_level_one(dim)
{
break;
}
}
level += 1;
}
iterator.reset_to_left_level_zero(dim);
Ok(())
}