Multigrid extension for NeuroFEM 
  last change: 04.03.2008  

New solver - Initialization methods

The initialization of the FMG solver is done inside the Init method. This method gets some parameters from NeuroFEM which describe the problem and the locations where informations are stored (right hand side and solution of the problem). The job of the initalization method is to store the parameters, to acquire informations about the problem, to build the grid and to set the local stiffness matrices.
Code from file "neurofem/fmg/fmgsolver.cpp"
void FMGSolver::Init(int ElemPerCell_, int numUnknowns_, int numCells_, doublereal* rhs_, doublereal* sol_) {
  // Store the variables
  ElemPerCell = ElemPerCell_;
  numUnknowns = numUnknowns_;
  numCells = numCells_;
  rhs = rhs_;
  sol = sol_;
  // Calculate the bounding box
  GetBoundingBox();
  // Calculate padding
  sizeX = PaddingPower2(NumCellsX);
  sizeY = PaddingPower2(NumCellsY);
  sizeZ = PaddingPower2(NumCellsZ);
  // Building up the grids etc.
  BuildGrids();
  // Assign the local stiffness matrices
  BuildMatrix();
}

Part of the initalization method is the BuildGrids method. The task of this method is to create ParExPDE objects describing the problem to be solved. This includes building up a hierarchical grid structure, allocating some variables and do some work for normalizing residual.
Code from file "neurofem/fmg/fmgsolver.cpp"
void FMGSolver::BuildGrids() {
  int px = sizeX;
  int py = sizeY;
  int pz = sizeZ;
  int levelx, levely, levelz, help;
  // Find the maximal level
  help = px; levelx = 1;
  while (help >= 8) {
    levelx++;
    help /= 2;
  }
  help = py; levely = 1;
  while (help >= 8) {
    levely++;
    help /= 2;
  }
  help = pz; levelz = 1;
  while (help >= 8) {
    levelz++;
    help /= 2;
  }
  number_levels = levelx;
  if (levely > number_levels) number_levels = levely;
  if (levelz > number_levels) number_levels = levelz;
  partition = ParExPDE::PartitionFactory::createMultigridCuboid(1.0, 1.0, 1.0, px+1, py+1, pz+1, number_levels);
  // *** Building some variables ***
  transfer_operator = new ParExPDE::StandardFETransferOperator;
  f = new ParExPDE::ETVariable(&(*partition), &(*transfer_operator));
  u = new ParExPDE::ETVariable(&(*partition), &(*transfer_operator));
  r = new ParExPDE::ETVariable(&(*partition), &(*transfer_operator));
  // *** Build the normalization ***
  normi = new double[number_levels];
  for ( int level = 0; level < number_levels - 1; level++ ) {
    *u = 1.0 | ParExPDE::interior_points;
    normi[ level ] = ParExPDE::ScalarProduct( *u, *u, &ParExPDE::interior_points );
    (*u).levelDown();
  }
  *u = 1.0 | ParExPDE::interior_points;
  normi[ number_levels - 1 ] = ParExPDE::ScalarProduct( *u, *u, &ParExPDE::interior_points );
  for ( int level = 0; level < number_levels - 1; level++ ) (*u).levelUp();
}

A second method used in the Init method is the BuildMatrix method. This method is used to build up the system matrix by transfering local stiffness matrices from NeuroFEM to ParExPDE. Additional this method needs to set the cells added as bounding box to an isolating state. The local stiffness matrices inside of NeuroFEM can be accessed by the volele_ method which requires a whole bunch of parameters. Since this method is very complicated, no detailed describtion of this method can be given here.
Code from file "neurofem/fmg/fmgsolver.cpp"
void FMGSolver::BuildMatrix() {
  ParExPDE::LocalStiffnessMatrix* lsm = op.getLocalStiffnessMatrix(&(*partition));
  ParExPDE::CuboidStorageBase* lsm_cells = lsm->getVariable()->getCuboidArray(0);
  // Calculate the isolator
  doublereal initls[8][8];
  for (int p = 0; p < 8; p++) for (int q = 0; q < 8; q++) {
    if (p == q) initls[p][q] = 1e-16;
    else initls[p][q] = 0.0;
  }
  // Initialize the localstiffness
  for (int k = 0; k < lsm_cells->getSize(2); k++) for (int j = 0; j < lsm_cells->getSize(1); j++) for (int i = 0; i < lsm_cells->getSize(0); i++) {
    ParExPDE::LocalStiffnessMatrixCell& cell = lsm_cells->valueAt(i,j,k);
    for (int p = 0; p < 8; p++) for (int q = 0; q < 8; q++) cell.valueAt(p,q) = initls[p][q];
  }
  // Set the real stiffness matrices
  int i,j,k;
  doublereal localstiff[ElemPerCell*ElemPerCell];
  doublereal dummy1[ElemPerCell][ElemPerCell];
  doublereal dummy2[ElemPerCell][ElemPerCell];
  for (int l = 1; l < numCells; l++) {
    MapCellToIndices(l, &i,&j,&k);
    ParExPDE::LocalStiffnessMatrixCell& cell = lsm_cells->valueAt(i,j,k);
    // Calculating local stiffness matrix
    volele_((integer *) VolumeGrid.getGridElementNodes(),
    (integer *) VolumeGrid.getGridElementType(),
    (integer *) VolumeGrid.getGridElementNumberNodes(),
    (doublereal *) VolumeGrid.getGridNodePosition(),
    (doublereal *) VolumeGrid.getGridElementConductivity(),
    (doublereal *) localstiff,
    (doublereal *) dummy1,
    (doublereal *) dummy2,
    (integer *) &numCells,
    (integer *) &numUnknowns,
    (integer *) &l);
    // Operator symmetric (no fortran reordering required)
    for (int p = 0; p < 8; p++) for (int q = 0; q < 8; q++) cell.valueAt(p,q) = localstiff[p*ElemPerCell+q];
  }
  op.setLSMInitialized(&(*partition));
}