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

Draft / notes / memo

By applying the harmonic stencil to itself (i.e. the Laplace operator):

\[ \Delta f(x,y) \approx \frac{f(x-h,y) + f(x+h,y) + f(x,y-h) + f(x,y+h) - 4f(x,y)}{h^2}, \]with the convolution operator we obtain the biharmonic stencil.

harmonic and biharmonic stencil (local host)

or here

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

Biharmonic distance (based on cotan weights)

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

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

Here we apply the biharmonic stencil iteratively, it should converge to biharmonic weights.

template< typename T> void diffuse_biharmonic(Grid2_ref grid, Grid2_const_ref mask, int nb_iter = 128*16) { tbx_assert( grid.size() == mask.size()); Vec2_i neighs_1[4] = {Vec2_i( 1, 0), Vec2_i( 0, 1), Vec2_i(-1, 0), Vec2_i( 0,-1)}; Vec2_i neighs_2[4] = {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 // to 'mask' int mask_nb_elts = mask.size().product(); std::vector idx_list; idx_list.reserve( mask_nb_elts ); 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; } } }

No comments