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

harmonic and biharmonic stencil (local host)

or here

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

Biharmonic distance (based on cotan weights)

Biharmonic Coordinates

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.

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

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: