# 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.

Note: finding the harmonic stencil can be simply done by summing the 2nd derivative central finite difference (first order accuracy) Advanced example of application of the harmonic stencil to solve a PDE or here

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

Biharmonic distance (based on cotan weights)

To look into as well: FinDiff is a python package that computes finite difference stencils and coefficients in any dimension, and can numerically compute derivatives or solve PDEs.

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,
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;
}
}
}



Related: 3D laplacian Just a comment: In Finite Differences, it can be much convenient to write it in terms of Kronecker products: if $A$ is the classical $[1,-2,1]$ matrix, you have that the 3D Laplacian is $$A = I \otimes I \otimes A + I \otimes A \otimes I + A \otimes I \otimes I$$ where $I$ is the identity matrix.

Paper to check out: Automated and Parallel Code Generation for Finite-Differencing Stencils with Arbitrary Data Types which describes the 2d and 3d laplacian and bilaplacian.

One comment

Hi!
What book is the “harmonic and biharmonic stencil” link from?
——————————-
Rodolphe: “Numerical Solution of the Biharmonic Equation” Petter E. Borstad (Phd thesis)

Vassillen Chizhov - 03-06-’22 20:55
(optional field, I won't disclose or spam but it's necessary to notify you if I respond to your comment)
All html tags except <b> and <i> will be removed from your comment. You can make links by just typing the url or mail-address.
Anti-spam question: 