36,38d35
< static PetscErrorCode MyInterpolation(Mat mat, Vec coarse, Vec fine);
< static PetscErrorCode MyRestriction(Mat mat, Vec fine, Vec coarse);
< 
52d48
<   PetscCall(PetscPrintf(PETSC_COMM_WORLD,">> Created DMshell %p (%p)\n",*shell,da));
77d72
<   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Calling KSPSolve from main\n"));
105c100
<   DM da1, da2;    // dm1 è la coarse, dm2 è la fine
---
>   DM da1, da2;
110,157c105
<   // PetscCall(DMCreateInterpolation(da1, da2, mat, vec));
<   PetscCall(PetscPrintf(PETSC_COMM_WORLD,">> Create interpolation from %p(%p) to %p(%p)\n",dm1,da1,dm2,da2));
< 
<   // getting the corners of the grids
<   PetscInt xLenCoarse, xLenFine;  
<   PetscCall(DMDAGetCorners(da1, NULL, NULL, NULL, &xLenCoarse, NULL, NULL));
<   PetscCall(DMDAGetCorners(da2, NULL, NULL, NULL, &xLenFine, NULL, NULL));
<   
<   PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, xLenFine, xLenCoarse, NULL, mat));
<   PetscCall(MatShellSetOperation(*mat, Mstatic ATOP_MULT, (void(*)(void))MyInterpolation));
<   PetscCall(MatShellSetOperation(*mat, MATOP_MULT_TRANSPOSE, (void(*)(void))MyRestriction));
<   PetscCall(MatSetUp(*mat));
< 
< 
<   PetscCall(MatViewFromOptions(*mat, NULL, "-int_view"));
<   //PetscCall(VecViewFromOptions(*vec, NULL, "-int_view"));
<   PetscFunctionReturn(PETSC_SUCCESS);
< }
< 
< static PetscErrorCode MyInterpolation(Mat mat, Vec coarse, Vec fine){
<   PetscFunctionBeginUser;
< 
<   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Inside myInterpolation\n"));
< 
<   PetscInt size;
< 
<   // get size of the vector
<   PetscCall(VecGetSize(coarse, &size));
< 
<   PetscCall(VecViewFromOptions(coarse, NULL, "-int_view"));
<   PetscCall(VecViewFromOptions(fine, NULL, "-int_view"));
< 
<   // setup the two vectors
<   const PetscScalar *vecCoarse;
<   PetscScalar *vecFine;
<   PetscCall(VecGetArrayRead(coarse, &vecCoarse));
<   PetscCall(VecGetArray(fine, &vecFine));
< 
<   for (PetscInt i = 0; i < size; i++)
<   {
<     if (i%2 == 0) vecFine[i] = vecCoarse[i];
< 
<     else vecFine[i] = 0.5*(vecCoarse[i-1]+vecCoarse[i+1]);
<   }
<   
<   PetscCall(VecRestoreArrayRead(coarse, &vecCoarse)); 
<   PetscCall(VecRestoreArray(fine, &vecFine));
< 
---
>   PetscCall(DMCreateInterpolation(da1, da2, mat, vec));
164,166c112
<   // Mat tmat;
<   
<   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"<< Create restriction\n"));
---
>   Mat tmat;
171,214c117,119
<   // PetscCall(DMCreateInterpolation(da1, da2, &tmat, NULL));
<   // PetscCall(MatTranspose(tmat, MAT_INITIAL_MATRIX, mat));
<   // PetscCall(MatDestroy(&tmat));
< 
<   // getting the corners of the grids
<   PetscInt xLenCoarse, xLenFine;  
<   PetscCall(DMDAGetCorners(da1, NULL, NULL, NULL, &xLenCoarse, NULL, NULL));
<   PetscCall(DMDAGetCorners(da2, NULL, NULL, NULL, &xLenFine, NULL, NULL));
<     
<   PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, xLenCoarse, xLenFine, NULL, mat));
<   PetscCall(MatShellSetOperation(*mat, MATOP_MULT, (void(*)(void))MyRestriction));
<   PetscCall(MatShellSetOperation(*mat, MATOP_MULT_TRANSPOSE, (void(*)(void))MyInterpolation));
< 
<   PetscCall(MatViewFromOptions(*mat, NULL, "-int_view")); 
<   PetscFunctionReturn(PETSC_SUCCESS);
< }
< 
< static PetscErrorCode MyRestriction(Mat mat, Vec fine, Vec coarse){
<   PetscFunctionBeginUser;
< 
<   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Inside myRestriction\n"));
< 
<   PetscInt size;
< 
<   // get size of the vector
<   PetscCall(VecGetSize(coarse, &size));
< 
<   PetscCall(VecViewFromOptions(coarse, NULL, "-int_view"));
<   PetscCall(VecViewFromOptions(fine, NULL, "-int_view"));
< 
<   // setup the two vectors
<   const PetscScalar *vecFine;
<   PetscScalar *vecCoarse;
<   PetscCall(VecGetArray(coarse, &vecCoarse));
<   PetscCall(VecGetArrayRead(fine, &vecFine));
< 
<   for (PetscInt i = 0; i < size; i++)
<   {
<     vecCoarse[i] = vecFine[i];
<   }
<   
<   PetscCall(VecRestoreArray(coarse, &vecCoarse)); 
<   PetscCall(VecRestoreArrayRead(fine, &vecFine));
< 
---
>   PetscCall(DMCreateInterpolation(da1, da2, &tmat, NULL));
>   PetscCall(MatTranspose(tmat, MAT_INITIAL_MATRIX, mat));
>   PetscCall(MatDestroy(&tmat));
245d149
<   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Inside Refine\n"));
257d160
<   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Inside Coarsen\n"));
274d176
<   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "computeRHS on grid %d\n",mx));
297d198
<   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "computeMatrix on grid %d\n",mx));