[petsc-users] Howto evaluate function on grid

[Klaus Zimmermann]

Dear all,

I want to evaluate a function on a grid. Right now I have some custom 
code (see the end of this message).
I am now thinking of rewriting this using more PETSC facilities because 
in the future we want to extend the
program to higher dimensions and also unstructured grids. As I 
understand it now there are several
candidates in PETSC for doing this:

1) PF
2) DALocalFunction
3) FIAT (?)

Could you please advise on what should be used?
An additional problem is that several distinct functions should be 
evaluated at the same time due to the
due to the reuse of intermediate results.

Any help is appreciated!
Thanks in advance,
Klaus

------------------8<------------------8<------------------8<------------------8<------------------8<--------------

PetscErrorCode EvaluatePsiAndGradPsi(AppCtx *user) {
   PetscErrorCode ierr;
   PetscInt i, j, cxs, cxm, cys, cym;
   PetscScalar **lpsi;
   PetscScalar **lGradPsi_x, **lGradPsi_y;
   Vec gc;
   DACoor2d **coors;
   DA cda;
   ierr = DAGetCoordinateDA(user->zGrid, &cda);CHKERRQ(ierr);
   ierr = DAGetCoordinates(user->zGrid, &gc);CHKERRQ(ierr);
   ierr = DAVecGetArray(cda, gc, &coors);CHKERRQ(ierr);
   ierr = DAGetCorners(cda, &cxs, &cys, PETSC_NULL, &cxm, &cym, 
PETSC_NULL);CHKERRQ(ierr);

   ierr = DAVecGetArray(user->zGrid, user->psi, &lpsi);CHKERRQ(ierr);
   ierr = DAVecGetArray(user->zGrid, user->gradPsi_x, 
&lGradPsi_x);CHKERRQ(ierr);
   ierr = DAVecGetArray(user->zGrid, user->gradPsi_y, 
&lGradPsi_y);CHKERRQ(ierr);
   for(i=cys; i<cys+cym; i++) {
     for(j=cxs; j<cxs+cxm; j++) {
       PetscReal x = PetscRealPart(coors[i][j].x-coors[i][j].y),
     y = PetscRealPart(coors[i][j].y);
       if(x>=0) {
     ierr = EvaluatePsiAndGradPsi(user, x, y,
&(lpsi[i][j]),
&(lGradPsi_x[i][j]),
&(lGradPsi_y[i][j]));CHKERRQ(ierr);
       }
     }
   }
   ierr = DAVecRestoreArray(user->zGrid, user->gradPsi_x, 
&lGradPsi_x);CHKERRQ(ierr);
   ierr = DAVecRestoreArray(user->zGrid, user->gradPsi_y, 
&lGradPsi_y);CHKERRQ(ierr);
   ierr = DAVecRestoreArray(user->zGrid, user->psi, &lpsi);CHKERRQ(ierr);
   ierr = DAVecRestoreArray(cda, gc, &coors);CHKERRQ(ierr);
   ierr = VecDestroy(gc);CHKERRQ(ierr);
   ierr = DADestroy(cda);CHKERRQ(ierr);
   return 0;
}