edu.jhu.cs.cisst.algorithms.optimize.fmg.pde
Interface PDE

All Known Implementing Classes:

CachingPDE, VectorFieldDiffusion2D

public interface PDE

The interface to a linear elliptic partial differential equation (PDE) as needed by the Full Multigrid (FMG) algorithm and relaxation algorithms used by the FMG algorithm.

In the framework of finite differencing, every very linear elliptic PDE can be written as a system of linear equations in matrix form, namely Au = f, where u is the vector of the solution values, f is the vector of the right hand side (RHS) values (often termed the source term), and A is the matrix that encapsulates all the details of the analytic form of the PDE, finite differencing and grid size.

A one-dimensional example may illustrate this: The analytical PDE is (d^2 u(x))/(d x^2) = f(x). Finite-differencing is done on the regular grid x[j] = x[0] + i*h, where for instance at level 2 the size of the grid would be 2^2 + 1 = 5. The interior grid elements are identified by i = 1,2,3, whereas the boundary grid elements are identified by i = 0,4. The matrix form of the PDE will be written for the total range 0 <= i <= 4. h is the grid spacing. The discretized PDE becomes (u[i - 1] - 2*u[i] + u[i + 1])/(h^2) = f[i]. So the RHS vector f would be (f[0],f[1],f[2],f[3],f[4]), but the boundary values f[0] and f[4] are not needed for solving the PDE and are just included here so that the size of the RHS vector is correct. The solution vector u would be (u[0],u[1],u[2],u[3],u[4]), but the boundary values u[0] and u[4] are given by the boundary conditions that are used to solve the PDE. In the case of fixed (Dirichlet) boundary conditions, the values u[0] and u[4] are given before the PDE can be solved and are thus not a result of the solving process. In the case of periodic boundary conditions u[0] is at all times set equal to u[3] and u[4] is at all times set equal to u[1] so that u[0] and u[4] are indirectly a result of the solving process. The matrix A written as a vector of row vectors would then be 1/(h^2)*((-,-,-,-,-),(1,-2,1,0,0),(0,1,-2,1,0),(0,0,1,-2,1),(-,-,-,-,-)). The first and last row of A are meaningless because the matrix equation Au = f is never directly solved for the first and last value of u (i.e. u[0] and u[4]). Furthermore, the left hand side (u[i - 1] - 2*u[i] + u[i + 1])/(h^2) of the discretized PDE and the knowledge of the value of h is a much more compact representation of the matrix A than the matrix representation given above.

Higher-dimensional PDEs are treated equivalently. This can be achieved by defining a linear order of the grid elements. E.g., the order ((0,0),(0,1),(0,2),(1,0),(1,1),(1,2),(2,0),(2,1),(2,2)) of a 2-dimensional grid of size 3x3 leads to the RHS vector f (f[0][0],f[0][1],f[0][2],f[1][0],f[1][1],f[1][2],f[2][0],f[2][1],f[2][2]).

This interface uses some of the terminology explained above but does not include the notion of a linear order of grid element. Insteads, the methods defined herein always operate directly on 3-dimensional grids.

Author:

Gerald Loeffler (Gerald.Loeffler@univie.ac.at)

Method Summary

double	evaluate(ConstBoundaryGrid u, ConstNoBoundaryGrid f, int x, int y, int z) evaluate the solution u of the PDE for a given interior grid element on a grid of a specific size.
double	evaluateLHS(ConstBoundaryGrid u, ConstNoBoundaryGrid f, int x, int y, int z) evaluate the left hand side of the discretized PDE for a given interior grid element on a grid of a specific size.
double	getGridSpacing(int level) return the grid spacing for a grid of the given size.
NoBoundaryGrid	sampleRHS(int sx, int sy, int sz, int level) sample the right hand side f on a cubic grid of the given size.
BoundaryGrid	solve(ConstBoundaryGrid u, ConstNoBoundaryGrid f) Solve PDE directly at coarsest grid level.

Method Detail

sampleRHS

NoBoundaryGrid sampleRHS(int sx,
                         int sy,
                         int sz,
                         int level)

sample the right hand side f on a cubic grid of the given size.

It is not necessary to know the value of f at the boundary of the grid, so the return type of this method is of type NoBoundaryGrid. But note that the size of the grid passed to this method includes (as always) the boundary. E.g., a 65x65x65 grid would be described by a size of 65 and its interior values (the one to be set by this method) would comprise the grid elements (1,1,1) to (63,63,63). The boundary elements (where at least one index is either 0 or 64) are non-existent in a NoBoundaryGrid.

Parameters:

sx - the sx

sy - the sy

sz - the sz

level - the level

Returns:

the right hand side (source term) f sampled on a cubic grid of the given size

getGridSpacing

double getGridSpacing(int level)

return the grid spacing for a grid of the given size.

The grid spacing h gives the distance of two neighbouring grid elements in one dimension measured in real world distance units. It can thus only be calculated if the mapping from grid space to the real world is known. The implementer of this method is supposed to have this knowledge.

Parameters:

level - the level

Returns:

the grid spacing for a grid of the given size and a specific mapping from grid space to the real world

evaluateLHS

double evaluateLHS(ConstBoundaryGrid u,
                   ConstNoBoundaryGrid f,
                   int x,
                   int y,
                   int z)

evaluate the left hand side of the discretized PDE for a given interior grid element on a grid of a specific size.

The left hand side Au of the PDE is to be evaluated. For the 1-dimensional example discussed above this would be the expression (u[i - 1] - 2*u[i] + u[i + 1])/(h^2) for a given value of i. In our (3-dimensional) case the grid element is of course identified by 3 integer indices.

Parameters:

u - the entire (approximate) solution values sampled on a grid of a specific size

x - ,y,z the integer indices of the interior grid element at which the left hand side is to be evaluated (0 < x,y,z < (grid size - 1))

f - the f

y - the y

z - the z

Returns:

the value of the left hand side of the PDE at the specified position

evaluate

double evaluate(ConstBoundaryGrid u,
                ConstNoBoundaryGrid f,
                int x,
                int y,
                int z)

evaluate the solution u of the PDE for a given interior grid element on a grid of a specific size.

The discretized form of the PDE can analytically be solved for a grid element at a specific position. For the 1-dimensional example discussed above this would be the expression 1/2*(u[i - 1] + u[i + 1] - f[i]*(h^2)) for a given value of i. In our (3-dimensional) case the grid element is of course identified by 3 integer indices.

Parameters:

u - the entire (approximate) solution values sampled on a grid of a specific size

f - the entire right hand side (source term) f sampled on a grid of the same size

x - ,y,z the integer indices of the interior grid element at which the left hand side is to be evaluated (0 < x,y,z < (grid size - 1))

y - the y

z - the z

Returns:

the value of the left hand side of the PDE at the specified position

solve

BoundaryGrid solve(ConstBoundaryGrid u,
                   ConstNoBoundaryGrid f)

Solve PDE directly at coarsest grid level.

Parameters:

u - the u

f - the f

Returns:

the boundary grid