Actual source code: ex6.c


  2: static char help[] = "Creates a matrix using 9 pt stencil, and uses it to test MatIncreaseOverlap (needed for additive Schwarz preconditioner). \n\
  3:   -m <size>       : problem size\n\
  4:   -x1, -x2 <size> : no of subdomains in x and y directions\n\n";

  6: #include <petscksp.h>

  8: PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
  9: {
 10:   Ke[0]  = H / 6.0;
 11:   Ke[1]  = -.125 * H;
 12:   Ke[2]  = H / 12.0;
 13:   Ke[3]  = -.125 * H;
 14:   Ke[4]  = -.125 * H;
 15:   Ke[5]  = H / 6.0;
 16:   Ke[6]  = -.125 * H;
 17:   Ke[7]  = H / 12.0;
 18:   Ke[8]  = H / 12.0;
 19:   Ke[9]  = -.125 * H;
 20:   Ke[10] = H / 6.0;
 21:   Ke[11] = -.125 * H;
 22:   Ke[12] = -.125 * H;
 23:   Ke[13] = H / 12.0;
 24:   Ke[14] = -.125 * H;
 25:   Ke[15] = H / 6.0;
 26:   return 0;
 27: }
 28: PetscErrorCode FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
 29: {
 30:   r[0] = 0.;
 31:   r[1] = 0.;
 32:   r[2] = 0.;
 33:   r[3] = 0.0;
 34:   return 0;
 35: }

 37: int main(int argc, char **args)
 38: {
 39:   Mat         C;
 40:   PetscInt    i, m = 2, N, M, idx[4], Nsub1, Nsub2, ol = 1, x1, x2;
 41:   PetscScalar Ke[16];
 42:   PetscReal   h;
 43:   IS         *is1, *is2, *islocal1, *islocal2;
 44:   PetscBool   flg;

 47:   PetscInitialize(&argc, &args, (char *)0, help);
 48:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 49:   N  = (m + 1) * (m + 1); /* dimension of matrix */
 50:   M  = m * m;             /* number of elements */
 51:   h  = 1.0 / m;           /* mesh width */
 52:   x1 = (m + 1) / 2;
 53:   x2 = x1;
 54:   PetscOptionsGetInt(NULL, NULL, "-x1", &x1, NULL);
 55:   PetscOptionsGetInt(NULL, NULL, "-x2", &x2, NULL);
 56:   /* create stiffness matrix */
 57:   MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 9, NULL, &C);

 59:   /* forms the element stiffness for the Laplacian */
 60:   FormElementStiffness(h * h, Ke);
 61:   for (i = 0; i < M; i++) {
 62:     /* node numbers for the four corners of element */
 63:     idx[0] = (m + 1) * (i / m) + (i % m);
 64:     idx[1] = idx[0] + 1;
 65:     idx[2] = idx[1] + m + 1;
 66:     idx[3] = idx[2] - 1;
 67:     MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES);
 68:   }
 69:   MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
 70:   MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);

 72:   for (ol = 0; ol < m + 2; ++ol) {
 73:     PCASMCreateSubdomains2D(m + 1, m + 1, x1, x2, 1, 0, &Nsub1, &is1, &islocal1);
 74:     MatIncreaseOverlap(C, Nsub1, is1, ol);
 75:     PCASMCreateSubdomains2D(m + 1, m + 1, x1, x2, 1, ol, &Nsub2, &is2, &islocal2);

 77:     PetscPrintf(PETSC_COMM_SELF, "flg == 1 => both index sets are same\n");
 78:     if (Nsub1 != Nsub2) PetscPrintf(PETSC_COMM_SELF, "Error: No of index sets don't match\n");

 80:     for (i = 0; i < Nsub1; ++i) {
 81:       ISEqual(is1[i], is2[i], &flg);
 82:       PetscPrintf(PETSC_COMM_SELF, "i =  %" PetscInt_FMT ",flg = %d \n", i, (int)flg);
 83:     }
 84:     for (i = 0; i < Nsub1; ++i) ISDestroy(&is1[i]);
 85:     for (i = 0; i < Nsub2; ++i) ISDestroy(&is2[i]);
 86:     for (i = 0; i < Nsub1; ++i) ISDestroy(&islocal1[i]);
 87:     for (i = 0; i < Nsub2; ++i) ISDestroy(&islocal2[i]);

 89:     PetscFree(is1);
 90:     PetscFree(is2);
 91:     PetscFree(islocal1);
 92:     PetscFree(islocal2);
 93:   }
 94:   MatDestroy(&C);
 95:   PetscFinalize();
 96:   return 0;
 97: }

 99: /*TEST

101:    test:
102:       args: -m 7

104: TEST*/