Question about divergence stencils

[ghost]

Hello,

What are the difference between the functions GetDivergence1 and GetDivergence2 in MultiGridOctreeData.System.inl?

GetDivergence1 and GetDivergence2 refer to the inner product tables dv and vd, respectfully. These tables only differ in sign, due to the B-Spline component derivatives.

Some explanation is written along SetCentralDivergenceStencil:

// if( scatter ) normals come from the center node
// else          normals come from the neighbors
template< int Degree1 , int Degree2 >
void SystemCoefficients< Degree1 , Degree2 >::SetCentralDivergenceStencil( const typename FunctionIntegrator::Integrator& integrator , Stencil< Point3D< double > , OverlapSize >& stencil , bool scatter )
{
    int center = ( 1<<integrator.depth() )>>1;
    int offset[] = { center , center , center };
    for( int x=0 ; x<OverlapSize ; x++ ) for( int y=0 ; y<OverlapSize ; y++ ) for( int z=0 ; z<OverlapSize ; z++ )
    {
        int _offset[] = { x+center-OverlapEnd , y+center-OverlapEnd , z+center-OverlapEnd };
        stencil.values[x][y][z] = scatter ? GetDivergence1( integrator , _offset , offset ) : GetDivergence2( integrator , _offset , offset );
    }
}

Could you provide some further intuition behind these functions, and what scatter means?

[mkazhdan]

It's been a while. But if I recall, the difference is whether you are integrating the dot-product if the gradients of two functions or the product of one with the Laplacian of the other. By Stokes' Theorem these are the same, but the first only requires the functions to be weakly differentiable, so it works for first order B-splines. (The latter wouldn't work as cleanly because the Laplacian explodes and so you would have to track what goes on the boundary of the integration domain.)

[ghost]

There are 4 tables for the inner products of the B-spline components, with different combinations of offset. dd as the inner product <B, B> dv as <grad(B), B> vd as <B, grad(B)> vv as <grad(B), grad(B)>

I think my confusion is that the underlying tables dv and vd are the same thing. However, their signs may differ.

[mkazhdan]

It could be that for the second and third case, one of the B's is a vector field and the other is scalar field. In one case you are integrating the dot product of a vector field with the gradient of a scalar function. In the other you are integrating the product of the divergence of a vector field with a scalar function. Moving the differentiation from one side to the other should flip the sign.

[ghost]

In the snippet below (from Version 7.0), both tables <grad(B), B> and <B, grad(B)>, are associated with integrating the dot product of a vector field with the gradient of a scalar function.

When the GRADIENT_DOMAIN_SOLUTION flag is enabled, they are both associated with integrating the dot product of a vector field with the gradient of a scalar function, where GRADIENT_DOMAIN_SOLUTION formulates the constraints as <B, div(V)>, instead of <grad(B), V>.

template< class Real >
Point3D< double > Octree< Real >::GetDivergence1( const typename BSplineData< 2 >::Integrator& integrator , int d , const int off1[] , const int off2[] , bool childParent ) const
{
    double vv[] =
    {
        integrator.dot( d , off1[0] , off2[0] , false , false , childParent ) ,
        integrator.dot( d , off1[1] , off2[1] , false , false , childParent ) ,
        integrator.dot( d , off1[2] , off2[2] , false , false , childParent )
    };
#if GRADIENT_DOMAIN_SOLUTION
    // Take the dot-product of the vector-field with the gradient of the basis function
    double vd[] = 
    {
        integrator.dot( d , off1[0] , off2[0] , false , true , childParent ) ,
        integrator.dot( d , off1[1] , off2[1] , false , true , childParent ) ,
        integrator.dot( d , off1[2] , off2[2] , false , true , childParent )
    };
    return  Point3D< double >( vd[0]*vv[1]*vv[2] , vv[0]*vd[1]*vv[2] , vv[0]*vv[1]*vd[2] );
#else // !GRADIENT_DOMAIN_SOLUTION
    // Take the dot-product of the divergence of the vector-field with the basis function
    double dv[] = 
    {
        integrator.dot( d , off1[0] , off2[0] , true , false , childParent ) ,
        integrator.dot( d , off1[1] , off2[1] , true , false , childParent ) ,
        integrator.dot( d , off1[2] , off2[2] , true , false , childParent )
    };
    return  -Point3D< double >( dv[0]*vv[1]*vv[2] , vv[0]*dv[1]*vv[2] , vv[0]*vv[1]*dv[2] );
#endif // GRADIENT_DOMAIN_SOLUTION
}

template< class Real > 
Point3D< double > Octree< Real >::GetDivergence2( const typename BSplineData< 2 >::Integrator& integrator , int d , const int off1[] , const int off2[] , bool childParent ) const
{
    double vv[] =
    {
        integrator.dot( d , off1[0] , off2[0] , false , false , childParent ) ,
        integrator.dot( d , off1[1] , off2[1] , false , false , childParent ) ,
        integrator.dot( d , off1[2] , off2[2] , false , false , childParent )
    };
#if GRADIENT_DOMAIN_SOLUTION
    // Take the dot-product of the vector-field with the gradient of the basis function
    double dv[] = 
    {
        integrator.dot( d , off1[0] , off2[0] , true , false , childParent ) ,
        integrator.dot( d , off1[1] , off2[1] , true , false , childParent ) ,
        integrator.dot( d , off1[2] , off2[2] , true , false , childParent )
    };
    return  Point3D< double >( dv[0]*vv[1]*vv[2] , vv[0]*dv[1]*vv[2] , vv[0]*vv[1]*dv[2] );
#else // !GRADIENT_DOMAIN_SOLUTION
    // Take the dot-product of the divergence of the vector-field with the basis function
    double vd[] = 
    {
        integrator.dot( d , off1[0] , off2[0] , false , true , childParent ) ,
        integrator.dot( d , off1[1] , off2[1] , false , true , childParent ) ,
        integrator.dot( d , off1[2] , off2[2] , false , true , childParent )
    };
    return -Point3D< double >( vd[0]*vv[1]*vv[2] , vv[0]*vd[1]*vv[2] , vv[0]*vv[1]*vd[2] );
#endif // GRADIENT_DOMAIN_SOLUTION
}

• integrator.dot( d , off1[0] , off2[0] , false , false , childParent ) corresponds to <B, B>
• integrator ..., true, false, ...) corresponds to <grad(B), B>, etc.