Actual source code: ex6.c

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

  5: #include <petscksp.h>

  7: static PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
  8: {
  9:   PetscFunctionBeginUser;
 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:   PetscFunctionReturn(PETSC_SUCCESS);
 27: }

 29: #if 0
 30: // unused
 31: static PetscErrorCode FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
 32: {
 33:   PetscFunctionBeginUser;
 34:   r[0] = 0.;
 35:   r[1] = 0.;
 36:   r[2] = 0.;
 37:   r[3] = 0.0;
 38:   PetscFunctionReturn(PETSC_SUCCESS);
 39: }
 40: #endif

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

 51:   PetscFunctionBeginUser;
 52:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 53:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 54:   N  = (m + 1) * (m + 1); /* dimension of matrix */
 55:   M  = m * m;             /* number of elements */
 56:   h  = 1.0 / m;           /* mesh width */
 57:   x1 = (m + 1) / 2;
 58:   x2 = x1;
 59:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-x1", &x1, NULL));
 60:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-x2", &x2, NULL));
 61:   /* create stiffness matrix */
 62:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 9, NULL, &C));

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

 77:   for (ol = 0; ol < m + 2; ++ol) {
 78:     PetscCall(PCASMCreateSubdomains2D(m + 1, m + 1, x1, x2, 1, 0, &Nsub1, &is1, &islocal1));
 79:     PetscCall(MatIncreaseOverlap(C, Nsub1, is1, ol));
 80:     PetscCall(PCASMCreateSubdomains2D(m + 1, m + 1, x1, x2, 1, ol, &Nsub2, &is2, &islocal2));

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

 85:     for (i = 0; i < Nsub1; ++i) {
 86:       PetscCall(ISEqual(is1[i], is2[i], &flg));
 87:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "i =  %" PetscInt_FMT ",flg = %d \n", i, (int)flg));
 88:     }
 89:     for (i = 0; i < Nsub1; ++i) PetscCall(ISDestroy(&is1[i]));
 90:     for (i = 0; i < Nsub2; ++i) PetscCall(ISDestroy(&is2[i]));
 91:     for (i = 0; i < Nsub1; ++i) PetscCall(ISDestroy(&islocal1[i]));
 92:     for (i = 0; i < Nsub2; ++i) PetscCall(ISDestroy(&islocal2[i]));

 94:     PetscCall(PetscFree(is1));
 95:     PetscCall(PetscFree(is2));
 96:     PetscCall(PetscFree(islocal1));
 97:     PetscCall(PetscFree(islocal2));
 98:   }
 99:   PetscCall(MatDestroy(&C));
100:   PetscCall(PetscFinalize());
101:   return 0;
102: }

104: /*TEST

106:    test:
107:       args: -m 7

109: TEST*/