C.2 An Example

In the following I provide a C++ code snippet that initializes a FieldVar. References to an FieldVar and a Mesh are passed in. The FieldVar is assumed to hold 3 fields e.g temperature, pressure and concentration, which are all cell-centered and initialized using analytical functions.

#include <cassert>
#include ``math.h''

void initializeFieldVar( FieldVar &fv, Mesh &mesh )
{

  double Pi = 3.1416 ;

  // Make sure we're working in 2D
  assert( mesh.getDimension() == 2 );

  // get deltax, deltay
  double dx = mesh.getDistances()[0] / mesh.getShape()[0];
  double dy = mesh.getDistances()[1] / mesh.getShape()[1];

  // make sure that the incoming fieldvariable has 3 fields
  assert( fv.getNVars() == 3 );

  // make sure it is cell-centered
  assert( fv.getMeshColl() == CENTERS ) ;

  /*
    Time to start looping over all the patches on this processor
    stored on this proc. Figure out how many patches exist.
   */
   int npatches = fv.getPatchCount();

   // loop over patches
   for( int ipatch = 0; ipatch < npatches; ipatch++ )
   {
     // get the spatial starting locations of the cell-centered arrays. 
     double startX = fv.getLowerCorner(ipatch)[0] * dx + 0.5*dx ;
     double startY = fv.getLowerCorner(ipatch)[1] * dy + 0.5*dy ;

     // get the dimensions of the patch
     int nx = fv.getShape(ipatch)[0] ;
     int ny = fv.getShape(ipatch)[1] ;

     /*
       get the current time and fetch the data pointer for the array
       which is supposed to hold the data for the current time
      */
     int itime = fv.getCurrentTime() ; 
     double *data = (double *) fv.get_data_raw(itime, ipatch)

     // OK, all done. time to initialize. loop over all fields in 
     // the FieldVar.
     for( int ifield = 0; ifield < fv.getNVars(); ifield++ )
     {
       // loop over all points in this patch
       for ( int j = 0; j < ny; j++ )
         for ( int i = 0; i < nx; i++ )
         {
           // figure out the x and y of this cell
           double x = startX + i*dx ;
           double y = startY + j*dy ;

           // figure out the stride for this cell and field;
           int index = i + j*nx + ifield*nx*ny ;
	   
           if ( ifield == 0 )         // Temperature field
              data[index] = cos(Pi*x) ;
           else if ( ifield == 1 )       // pressure field
              data[index] = sin(Pi*x) ;
           else                          // concentration
              data[index] = tanh(x) ;
         }  // End of loop over points

      } // end of loop over fields
       
   } // End of loop over patches

   // all done; return

   return ;
}