Embedding multigrid algorithms in the PZ environment

The PZ environment was modified and extended to include algorithms which will simplify the development of multigrid algorithms with the PZ environment. These changes should not affect existing software.

There is distinguishment between

• Changes made to the existing classes : these include modification to existing methods (but without modifying the underlying behaviour), changes made to the datastructure and addition of new methods to accomodate the multigrid method. Considering the multigrid method require the use of grids of different levels, the resulting code is now much closer to h-adaptivity.
  • Geometric mapping classes
  • Computational element classes
  • Matrix classes
• Development of new classes, particular for the implementation of the multigrid method. These include
  • A sparse block matrix to implement the transfer matrix
  • A grid generator which generates grids of several levels for triangular and quadrilateral elements
  • A transformation to map integration points on the fine element to the master element space of the coarse element

Changes made to existing classes

Geometric Mapping classes

• The geometric elements have been extended to allow the user to provide information on the hierarchy tree of the elements. As such methods like
  • void TGeoEl::SetFather(TGeoEl *father);
    // identifies the object father as the father of the current element
  • TGeoEl *TGeoEl::Father();// returns the father of the current element
  • void TGeoEl::SetSubElements(VoidPtrVec &subel);
    // identifies the elements pointed to by the elements of subel as the subelements of the current element
  • int TGeoEl::Level();
    // returns the level of the geometric element
  • virtual void TGeoEl::BuildTransform(TGeoEl *father, TTransform &t) = 0;
    // accumulates the transformation of the jacobian which maps the current
    // master element space into the space of the master element of the father
• The classes TGeoEl1d, TGeoElQ2d, TGeoElT2d have been extended to include the implementation of
  void BuildTransform(TGeoEl *father, TTransform &t);
• TGeoGrid::FindElement(--) which finds a boundary element connected to a node now only returns an element with zero level. This allows to call the TGeoGrid::GetBoundaryElements(---) for grids containing more than one level of elements. This method will return a vector of level zero elements between two nodes on the boundary.

Computational Element Classes

• The computational element classes have been modified to
  1. Unify standard methods of finite element computations. At this stage, CalcStiff, EvaluateError, ProjectFlux are implemented in the derived classes, causing duplication of code. The data structure of the computational element classes has been modified to enable to implement the above methods at the base TCompEl class. This will reduce the code and facilitate the incorporation of new (3D) elements.
  2. Implement a method to compute the correspondence of the shapefunctions of a small element within a large element. The matrix which describes this relation can be used as a multigrid transfer matrix. Using the modified data structure mentioned above, this method was implemented in the base TCompEl class. This implies that the projection is implemented for all elements of the PZ environment simultaneously.
  3. The computational elements now include the order of interpolation of the shapefunctions on their sides as a data. Previously the order of interpolation was computed when needed. The new approach has several advantages :
     • The element does not need to access its neighbours during element stiffness computations. This is particularly important for parallel implementations, where neighbouring elements may reside on different processors
     • During the computation of the projection of the shapefunctions mentioned above, the neighbouring information is not available
     • The implementation of low order elements with high order side interpolation is now available
  4. Several methods have been included to declare the element more precisely to the environment. All elements now declare the total number of shapefunctions, the number of shapefunctions on each side, the order of interpolation, its dimension (1 - 2 - 3).
  5. The declaration of the Shape() method has been made uniform for all elements. Its implementation is still particular in each sub element.
  6. The pointer to the integration rule is the same for all elements.
  7. The inconsistency which existed between the dimensions of the derivative of the shapefunctions in one dimension and two dimensions has been removed. The dimension of the derivative comes first and the number of shape functions comes second.
     There may still be places in the program where the indexes need to be inverted
• The computational grid has been extended to include a method which computes a transfer matrix which interpolates the solution of the coarse grid onto the fine grid or injects the residual of the fine grid onto the coarse grid.
• The implementation of the computation of the transfer matrix is not supported on grids which contain superelements and will probably not work on superelements also.

Matrix classes

The MultAdd and Multiply methods of the matrix classes (TMatrix and derived) have been extended to allow multiplication with a specified stride. This extension is needed because the interpolation matrix computed by the computational grid and elements specifies the interpolation between shapefunctions. However the vectors corresponding to solutions and right hand sides of systems of differential equations group their variables. This extension will also be particularly useful for multiplying a mass matrix with a solution vector in dynamical problems.

A bug which existed in the MultAdd method of most matrix classes has been fixed.

New Classes Developed

TTransfer : A sparse block matrix class

The transfer matrix between a fine grid and a coarse grid is typically very sparse. Moreover, in high order interpolation, the sparsity of the matrix is block oriented : the equations associated with a side of an element, with the interior of the element and with the vertices of the element form tightly coupled blocks of equations.

The TTransfer class implements a yale sparse matrix type storage form of matrices. This is one of the most efficient storage patterns known to the author and still allows to implement computationally efficient matrix vector multiplications. Due to its particular storage pattern, the user must supply the elements one row at a time, but the order of the rows may be arbitrary. There are no methods to delete elements, nor to add elements in an arbitrary manner. Arbitrary elements cannot be accessed, but the matrix vector multiplication (with stride != 1) has been implemented.

The public methods of the TTransfer class are :

• class TTransfer : public TMatrix
  // this class implements a rectangular sparse block matrix
  // it is assumed that the data is entered one row at a time
  // the matrix structure cannot be modified after being defined
• TTransfer();
• TTransfer(TBlock &row, TBlock &col);
  // the sparse matrix blocks are defined by row, col
• virtual void Print( char *name = NULL, ostream &out = cout , MatrixOutputFormat form = EFormatted);
• void SetBlocks(TBlock &row, TBlock &col);
  // this operation will reset the matrix to zero
  // with no rows defined
• int HasRowDefinition(int row);
  // returns 1 if the row is defined (i.e. has column entries)
• void AddBlockNumbers(int row, TIntVec &colnumbers);
  // will specify the sparsity pattern of row
• void SetBlockMatrix(int row, int col, TFMatrix &mat);
  // sets the row,col block equal to matrix mat
  // if row col was not specified by AddBlockNumbers, an error
  // will be issued and exit
• void MultAdd( TFMatrix &x, TFMatrix &y, TFMatrix &z,
  REAL alpha, REAL beta, int opt = 0, int stride = 1);
  // multiplies the transfer matrix and puts the result in z

TTransform : a simple transformation class

The TTransform class maps a point in space to another point through a linear transformation followed by a translation. The dimension of the transformation can be one-, two- or three dimensional and needs to be specified at construction. Although the transformation class includes regular full matrix objects (type TFMatrix), it does not use dynamic memory allocation. It is therefore suited for use as local objects. Two methods are provided to modify a transformation : one in which the user specifies the transformation matrices and another where two transformations are multiplied to form the combined transformation. A method is provided to apply the transformation to a vector and put the result in a second vector. The vectors may not be the same. The public interface of the transformation class is :

• class TTransform {
• TTransform(int dim);
• TTransform(TTransform &tr);
• TTransform &operator=(TTransform &t);
• void SetMatrix(TFMatrix &mult, TFMatrix &sum);
  // sets the transformation matrices
• TTransform Mult(TTransform &right);
  // multiply the transformation object to the right with right
• void Apply(DoubleVec &in, DoubleVec &out);
  // transforms the vector

TMultGenGrid : a class to generated a geometric grid with several levels

TMultGenGrid is an extension of the TGenGrid class developed by M. Santana, which allows to generate a grid with several uniform refinements in one step. Considering the TGeoGrid class holds the entire structure of the refined grid, the TMultGenGrid class still works on a single grid. The public interface of the TMultGenGrid class is :

• class TMultGenGrid : public TGenGrid
• TMultGenGrid (TIntVec &nx, DoubleVec &x0, DoubleVec &x1, int numlev);
  // initializes the grid object using
  // nx to specify the number of elements in either direction
  // x0 as the lower left corner of the domain
  // x1 as the upper right corner of the domain
  / numlev as the number of levels to be generated
• virtual ~TMultGenGrid();
• virtual void SetElementType(int type);
  // sets the element type to rectangular (type == 0)
  // or triangular (type == 1)
• virtual short Read(TGeoGrid &malha);
  // generates the grid structure by calling the appropriate methods
• virtual void SetBC(TGeoGrid &gr, int side, int bc);
  // sets the boundary condition along side of the rectangular grid to bc
• virtual void SetBC(TGeoGrid &gr, DoubleVec &start, DoubleVec &end, int bc);
  // sets the boundary condition between points start and end along the boundary to bc
• virtual void Print( char *name = NULL, ostream &out = cout);
  // prints the grid data