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

Name

Email (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:
√4 + 40 = ?

There is no character, howsoever good and fine, but it can be destroyed by
ridicule, howsoever poor and witless. Observe the ass, for instance: his
character is about perfect, he is the choicest spirit among all the humbler
animals, yet see what ridicule has brought him to. Instead of feeling
complimented when we are called an ass, we are left in doubt.
-- Mark Twain, "Pudd'nhead Wilson's Calendar"