TransferToPETSc.cpp

Go to the documentation of this file.

#ifdef with_petsc

#include <stdlib.h>
#include <vector>
#include <basic.h>
#include <arrayfunctions.h>
#include <simulation_object.h>
#include <bandedmatrix.h>
#include <petscinterface_struct.h>

#undef __FUNCT__
#define __FUNCT__ "TransferVecToPETSc"

int TransferVecToPETSc( const double* const u, 
                        Vec Y, 
                        void* ctxt, 
                        const int sim_idx,
                        const int offset   )
{
  PETScContext* context = (PETScContext*) ctxt;
  SimulationObject* sim = (SimulationObject*) context->simobj;
  double* Yarr;

  PetscFunctionBegin;
  VecGetArray(Y,&Yarr);
  std::vector<int> index(sim[sim_idx].solver.ndims,0);
  ArrayCopynD(  sim[sim_idx].solver.ndims,
                u,
                (Yarr+offset),
                sim[sim_idx].solver.dim_local,
                sim[sim_idx].solver.ghosts,
                0,
                index.data(),
                sim[sim_idx].solver.nvars );
  VecRestoreArray(Y,&Yarr);

  PetscFunctionReturn(0);
}

int TransferMatToPETSc( void *J,    
                        Mat   A,    
                        void *ctxt   )
{
  BandedMatrix    *M = (BandedMatrix*) J;
  PetscErrorCode  ierr     = 0;
  int             bs = M->BlockSize, nbands = M->nbands, bs2 = bs*bs;

  for (int i=0; i<M->nrows_local; i++) {
    int     colind[nbands];
    double  val[bs][bs*nbands];
    for (int n=0; n<nbands; n++) {
      colind[n] = M->ncol[nbands*i+n];
      for (int p=0; p<bs; p++) {
        for (int q = 0; q<bs; q++) {
          val[p][n*bs+q] = M->data[i*nbands*bs2+n*bs2+p*bs+q];
        }
      }
    }
    MatSetValuesBlocked(A,1,&M->nrow[i],M->nbands,&colind[0],&val[0][0],INSERT_VALUES); 
  }

  MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);

  return(0);
}

#endif