hodge.h

Hodge Decomposition

Hodge Decomposition

Let \Omega^k the space of the k-forms on a closed compact n-dimensional manifold, any k-form \omega can be expressed as \displaystyle \omega = d \alpha + \delta \beta + \gamma for some k-1 form \alpha, a k+1 form \beta and \gamma \in \Omega^k. Where in the exterior calculus (see Exterior Algebra)

d is the differential exterior derivative
\delta is the codifferential defined by \star d \star

and \star is the Hodge operator.

Given a k-form \omega there are unique \alpha, \beta,\gamma (Hodge theorem, see H. Flanders, Differential Forms with Applications to the Physical Sciences (Dover Books on Mathematics)) which can be found by solving

\delta \omega = \delta d \alpha
d \omega = d \delta \beta
\gamma = \omega - d \alpha - \delta \beta

In R^2 a vector X = (u,v) is associated to a 1-form \omega = u dx + v dy

we solve \delta \omega = \delta d \alpha

\delta \omega = \star d \star \omega \rightarrow \delta \omega_X = \star d ( -v dx + u dy)
\delta \omega = \star (-v_y dy \wedge dx + u_x dx \wedge dy)
\delta \omega = u_x + v_y

As \omega_X is a 1-form that implies \alpha \in \Omega^0 a function of real values, we found then

d \alpha = \alpha_x dx + \alpha_y \ dy
\delta d \alpha = \star d \star d \alpha = \alpha_{xx} + \alpha_{yy}

The first equation is then \displaystyle \alpha_{xx} + \alpha_{yy} = u_x + v_y a Poisson equation for \alpha.

Structure for the solver :

omega,
gradb,
rotc,
h

are mandatory vectors giving the Hodge decomposition. If the scalar alpha and beta exist they will return the scalar fields.

#include "poisson.h"

struct Hodge {
  // mandatory
  vector omega;
  vector gradb;
  vector rotc;
  vector h;
  // optional
  scalar alpha;
  scalar beta;
};

Call the Hodge decomposition

trace
mgstats hodge (struct Hodge p)
{

 vector omega =  p.omega;
 vector gradb = p.gradb;
 vector rotc = p.rotc;
 vector h = p.h;

  // Local
 scalar b[];
 scalar alpha[];
 scalar beta[];

Check the existence of alpha et beta

 if (p.alpha.i) {
   alpha = p.alpha;
 }

 if (p.beta.i) {
   beta = p.beta;
 }

paranoia

boundary  ((scalar *){omega});

Filling the rhs of the Poisson equation

  foreach()
    {
#if order == 2
      b[] = (omega.x[1,0] - omega.x[-1,0])/2./Delta + (omega.y[0,1] - omega.y[0,-1])/2./Delta;
#else
      b[] =   (omega.x[-2,0] - 8.*omega.x[-1,0] + 8.*omega.x[1,0]  -omega.x[2,0])/12./Delta;
      b[] +=  (omega.y[0,-2] - 8.*omega.y[0,-1] + 8.*omega.y[0,1]  -omega.y[0,2])/12./Delta;
    }
#endif

Solve \nabla^2 \alpha = b using the Poisson solver

 mgstats s = poisson (alpha, b,tolerance = 1e-6);

Now we apply the same procedure to d \omega = d \delta \beta

d\omega = \left( -\frac{\partial\,u}{\partial y} + \frac{\partial\,v}{\partial x} \right) \mathrm{d} x\wedge \mathrm{d} y

and

d \star d \star \beta = \left( \frac{\partial^2\,\beta}{\partial x ^ 2} + \frac{\partial^2\,\beta}{\partial y ^ 2} \right)d x\wedge d y

we have still a Poisson equation

\displaystyle \nabla^2 \beta = b = \left( -\frac{\partial\,u}{\partial y} + \frac{\partial\,v}{\partial x} \right)

  foreach()
    {
#if order == 2
      b[] = (omega.y[1,0] - omega.y[-1,0])/2./Delta - (omega.x[0,1] - omega.x[0,-1])/2./Delta;
#else
      b[] =  (omega.y[-2,0] - 8.*omega.y[-1,0] + 8.*omega.y[1,0]  -omega.y[2,0])/12./Delta;
      b[] -= (omega.x[0,-2] - 8.*omega.x[0,-1] + 8.*omega.x[0,1]  -omega.x[0,2])/12./Delta;
    }
#endif
  boundary ({beta,b});

Solve Poisson problem

 s = poisson (beta, b,tolerance = 1e-6);

Finally from the scalar fields \alpha and \beta$ we compute the vectors \displaystyle gradb = d \alpha \displaystyle rotc = \delta \beta and \gamma by subtraction.

Compute vectors

foreach()
    foreach_dimension()
#if order == 2
    gradb.x[] = (alpha[1,0] - alpha[-1,0]) / 2./Delta ;
#else
  gradb.x[] = (alpha[-2,0] - 8*alpha[-1,0]+ 8.* alpha[1,0] - alpha[2,0]) / 12./Delta ;
#endif

  
struct { double x, y; } a = {-1., 1.};
 
foreach()
    foreach_dimension()
#if order == 2
    rotc.x[] =  a.x * (beta[0,1] - beta[0,-1]) / 2./Delta ;
  #else
  rotc.x[] =  a.x * (beta[0,-2] - 8*beta[0,-1]+ 8.* beta[0,1] - beta[0,2]) / 12./Delta ;
#endif


  boundary ( (scalar *){gradb,rotc});

Compute h

  foreach()
    foreach_dimension()
    h.x[] = omega.x[] - gradb.x[] - rotc.x[];

  boundary  ((scalar *){h});

  
    return s;  // Fix : better return 

}

2D Test