# 2D biharmonic stencil a.k.a bilaplacian operator

Draft work in progress

By applying the harmonic stencil to itself with the convolution operator we obtain the biharmonic stencil. or here

For triangle mesh not sure where to look, perhaps here:

Biharmonic distance (based on cotan weights)

Biharmonic Coordinates

Localized bi-Laplacian Solver on a Triangle Mesh and Its Applications

Todo: 1d biharmonic stencil (1 -4 6 -4 1) (https://en.wikipedia.org/wiki/Finite_difference_coefficient) which you can infer by applying the harmonic stencil (1 -2 1) to itself

template< typename T>
void diffuse_biharmonic(Grid2_ref grid,
int nb_iter = 128*16)
{
Vec2_i neighs_1 = {Vec2_i( 1, 0),
Vec2_i( 0, 1),
Vec2_i(-1, 0),
Vec2_i( 0,-1)};

Vec2_i neighs_2 = {Vec2_i( -1, -1),
Vec2_i( -1,  1),
Vec2_i(  1,  1),
Vec2_i(  1, -1)};

// First we store every grid elements that need to be diffused according
std::vector idx_list;
for(Idx2 idx(mask.size(), 0); idx.is_in(); ++idx) {
if( mask( idx ) == 1 )
idx_list.push_back( idx.to_linear() );
}

// Diffuse values with a Gauss-seidel like scheme
for(int j = 0; j < nb_iter; ++j)
{
// Look up grid elements to be diffused
for(unsigned i = 0; i < idx_list.size(); ++i)
{
Idx2 idx(grid.size(), idx_list[i]);

T sum  = T(0);
T mean = T(0);

for(int i = 0; i < 4; ++i)
{
Idx2 idx_neigh = idx + neighs_1[i];
if( idx_neigh.is_in() /*check if inside grid*/ ){
T weight = (T)8;
sum = sum + grid( idx_neigh ) * weight;
mean += weight;
}
}

for(int i = 0; i < 4; ++i)
{
Idx2 idx_neigh = idx + neighs_2[i];
if( idx_neigh.is_in() /*check if inside grid*/ ){
T weight = (T)-2;
sum = sum + grid( idx_neigh ) * weight;
mean += weight;
}
}

for(int i = 0; i < 4; ++i)
{
Idx2 idx_neigh = idx + neighs_1[i] * 2;
if( idx_neigh.is_in() /*check if inside grid*/ ){
T weight = (T)-1;
sum = sum + grid( idx_neigh ) * weight;
mean += weight;
}
}

sum = sum / (T)( mean );
grid( idx ) = sum;
}
}
} 