1 static char help[] = "Tests VecView() contour plotting for 2d DMDAs.\n\n"; 2 3 /* 4 MATLAB must be installed to configure PETSc to have MATLAB engine. 5 Unless you have specific important reasons for using the MATLAB engine, we do not 6 recommend it. If you want to use MATLAB for visualization and maybe a little post processing 7 then you can use the socket viewer and send the data to MATLAB via that. 8 9 VecView() on DMDA vectors first puts the Vec elements into global natural ordering before printing (or plotting) 10 them. In 2d 5 by 2 DMDA this means the numbering is 11 12 5 6 7 8 9 13 0 1 2 3 4 14 15 Now the default split across 2 processors with the DM is (by rank) 16 17 0 0 0 1 1 18 0 0 0 1 1 19 20 So the global PETSc ordering is 21 22 3 4 5 8 9 23 0 1 2 6 7 24 25 Use the options 26 -da_grid_x <nx> - number of grid points in x direction, if M < 0 27 -da_grid_y <ny> - number of grid points in y direction, if N < 0 28 -da_processors_x <MX> number of processors in x directio 29 -da_processors_y <MY> number of processors in x direction 30 */ 31 32 #include <petscdm.h> 33 #include <petscdmda.h> 34 35 int main(int argc, char **argv) 36 { 37 PetscMPIInt rank; 38 PetscInt M = 10, N = 8; 39 PetscBool flg = PETSC_FALSE; 40 DM da; 41 PetscViewer viewer; 42 Vec local, global; 43 PetscScalar value; 44 DMBoundaryType bx = DM_BOUNDARY_NONE, by = DM_BOUNDARY_NONE; 45 DMDAStencilType stype = DMDA_STENCIL_BOX; 46 #if defined(PETSC_HAVE_MATLAB) 47 PetscViewer mviewer; 48 PetscMPIInt size; 49 #endif 50 51 PetscFunctionBeginUser; 52 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 53 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "", 300, 0, 300, 300, &viewer)); 54 #if defined(PETSC_HAVE_MATLAB) 55 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 56 if (size == 1) PetscCall(PetscViewerMatlabOpen(PETSC_COMM_WORLD, "tmp.mat", FILE_MODE_WRITE, &mviewer)); 57 #endif 58 59 PetscCall(PetscOptionsGetBool(NULL, NULL, "-star_stencil", &flg, NULL)); 60 if (flg) stype = DMDA_STENCIL_STAR; 61 62 /* Create distributed array and get vectors */ 63 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, bx, by, stype, M, N, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da)); 64 PetscCall(DMSetFromOptions(da)); 65 PetscCall(DMSetUp(da)); 66 PetscCall(DMCreateGlobalVector(da, &global)); 67 PetscCall(DMCreateLocalVector(da, &local)); 68 69 value = -3.0; 70 PetscCall(VecSet(global, value)); 71 PetscCall(DMGlobalToLocalBegin(da, global, INSERT_VALUES, local)); 72 PetscCall(DMGlobalToLocalEnd(da, global, INSERT_VALUES, local)); 73 74 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 75 value = rank + 1; 76 PetscCall(VecScale(local, value)); 77 PetscCall(DMLocalToGlobalBegin(da, local, ADD_VALUES, global)); 78 PetscCall(DMLocalToGlobalEnd(da, local, ADD_VALUES, global)); 79 80 flg = PETSC_FALSE; 81 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_global", &flg, NULL)); 82 if (flg) { /* view global vector in natural ordering */ 83 PetscCall(VecView(global, PETSC_VIEWER_STDOUT_WORLD)); 84 } 85 PetscCall(DMView(da, viewer)); 86 PetscCall(VecView(global, viewer)); 87 #if defined(PETSC_HAVE_MATLAB) 88 if (size == 1) { 89 PetscCall(DMView(da, mviewer)); 90 PetscCall(VecView(global, mviewer)); 91 } 92 #endif 93 94 /* Free memory */ 95 #if defined(PETSC_HAVE_MATLAB) 96 if (size == 1) PetscCall(PetscViewerDestroy(&mviewer)); 97 #endif 98 PetscCall(PetscViewerDestroy(&viewer)); 99 PetscCall(VecDestroy(&local)); 100 PetscCall(VecDestroy(&global)); 101 PetscCall(DMDestroy(&da)); 102 PetscCall(PetscFinalize()); 103 return 0; 104 } 105 106 /*TEST 107 108 test: 109 requires: x 110 nsize: 2 111 args: -nox 112 output_file: output/empty.out 113 114 TEST*/ 115