DGtal 1.4.0

The Poisson equation can be written as:
\[ \Delta \phi = f \]
where \(\Delta\) is the Laplace operator, \(\phi\) the solution of the problem and \(f\) the input of the problem. The Poisson equation can be viewed as the heat equation when steady state is reached. \(f\) is then the constant spatial heating power applied to the structure and \(\phi\) is the temperature reached in steady state (up to a spatial constant). Here we choose \(f\) to be a Dirac delta function \(\delta\) to simulate a punctual heating on the structure.
Depending on border conditions, \(\Delta\) may have some null eigenvalues. In order to make it solvable, we use a regularized version of the Laplace operator \(\Delta_\mathrm{reg}\), defined below
\[ \Delta_\mathrm{reg} = \Delta + \lambda I \]
where \(I\) is the identity operator and \(\lambda\) is the regularization scalar, chosen small with respect to \(\Delta\) eigenvalues. Linear regularized solution are the same solution as least square problem solution of the non regularized problem. \(I\) is generated using DiscreteExteriorCalculus.identity.
One can use the general Laplace operator definition.
\[ \Delta = \star d \star d = \delta d \]
However, for convenience DiscreteExteriorCalculus provides a method to get the Laplace operator directly.
In this example, we create an 1D linear structure embedded in a 2D space and we solve a non regularized Poisson equation on this structure.
Boundary conditions effect are illustrated for the two classical boundary conditions.
Boundary conditions are changed from Neumann to Dirichlet conditions by adding dangling 1cells at the beginning and at the end of the structure. Snippets are taken from testLinearStructure.cpp.
First, an empty calculus structure is created and, using simple for loops, it is filled with some 0cells and 1cells to from a linear structure. When filled manually with DiscreteExteriorCalculus.insertSCell, one can pass the primal and dual sizes of each cell which default to 1. Temperature nodes are associated primal 0cells. Note that to enforce Neumann boundary conditions, the linear structure has to end with 0cells. Moreover some 1cells are inserted as negative cells to match their orientation with the orientation of the linear structure.
The input heat vector \(f\), which is a primal 0form in the DEC formulation, is then created and is given values of a Dirac pulse shifted at the right position.
The primal non regularized Laplace operator \(\Delta\) is generated using DiscreteExteriorCalculus.laplace. The primal exterior derivative from 0form to 1form will be used to compute the gradient solution.
Now the problem is fully defined and there one thing left to do: solving it. The resolution is done by DiscreteExteriorCalculusSolver. This class takes the actual linear solver used as the second template parameter. Any class that validates the CLinearAlgebraSolver concept can be wrapped by DiscreteExteriorCalculusSolver. If DGtal was compiled with eigen support enabled as described in DEC linear solver, some solvers will be available in the EigenLinearAlgebraBackend. Here, we will use the EigenLinearAlgebraBackend::SolverSparseQR solver. Once created the solver is given the operator using DiscreteExteriorCalculusSolver.compute. The input dirac 0form is passed to DiscreteExteriorCalculusSolver.solve, which return the solution of the problem.
Since the dirac input is null everywhere except at a single point, this means that the second derivative of the solution is null everywhere except at the dirac position. An analytic form can be expressed as a continuous piecewise quadratic function. Numerical values of the solution fit analytic values with at least a relative precision of 1e5.
Dirichlet boundary condition fixes value of 0forms to zero along borders of the structure. Two dangling 1cells are added at each end of the Neumann structure to switch to Dirichlet boundary condition. Since those 1cells are not connected to a 0cell on one of their border, this will simulate the presence of zerovalued 0cell in those places through enforcing Dirichlet boundary conditions.
The input Dirac can be used as in the Neumann case since no 0cells has been added to the structure.
Laplace operator needs to be rebuild, even if the code doesn't change.
Solving the problem is achieved by using the same code as for the Neumann case. This time the analytic solution is piecewise linear and take a constant null value at the border of the structure, as expect from the Dirichlet boundary condition.
Numerical values of the solution fit analytic values with at least a relative precision of 1e5.
In this example, we create a 2D ring structure and we solve a regularized Poisson equation on this structure. Snippets are taken from exampleDiscreteExteriorCalculusSolve.cpp. First the structure is created from a digital set.
Then we compute the dual regularized Laplace operator with \(\lambda = 0.01 \).
Input kform \(y\) is a Dirac delta positioned on the dual 0cell at coordinates \((2,5)\). DiscreteExteriorCalculus.getIndex return the index of the dual 0form container associated with the dirac position cell.
The following illustration represents the dual of the input. Hence, since we have specified a dirac on "dual 0cell", its primal representation attaches information to (primal) 2cells.
We try to solve the problem using EigenLinearAlgebraBackend::SolverSimplicialLLT, but DiscreteExteriorCalculusSolver.isValid reports an error. Underlying linear algebra solver DiscreteExteriorCalculusSolver.solver.info reports a numerical_error.
Since the first solver failed, let's use another one: EigenLinearAlgebraBackend::SolverSimplicialLDLT for example. This time DiscreteExteriorCalculusSolver.isValid is true after computing problem solution and the solution dual 0form contains the solution of the problem.
In this example, we create an embedded 2D structure from a DigitalSurface and we solve a dual Poisson equation with Neumann boundaries condition. Snippets are taken from exampleDECSurface.cpp.
First we load the Alcapone image and extract its boundary surface. Note that the Khalimsky space domain opens the surface under the feet of Alcapone.
The DEC structure is then create using DiscreteExteriorCalculusFactory::createFromNSCells. Borders are not added to the structure.
Dual 0form \(\rho\) is the input of the Poisson problem. Left foot has positive value and right foot has negative value.
Dual Laplace operator is computed and passed to the solver using DiscreteExteriorCalculusSolver.compute. The dual 0form solution \(\phi\) is returned by DiscreteExteriorCalculusSolver.solve.