code/ex49.c

   1 //  /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 4    ./ex49 -mx 900 -my 900 -ksp_type fgmres
   2
   3
   4 static char help[] =  "   Solves the compressible plane strain elasticity equations in 2d on the unit domain using Q1 finite elements. \n\
   5    Material properties E (Youngs moduls) and nu (Poisson ratio) may vary as a function of space. \n\
   6    The model utilisse boundary conditions which produce compression in the x direction. \n\
   7 Options: \n\
   8      -mx : number elements in x-direciton \n\
   9      -my : number elements in y-direciton \n\
  10      -c_str : indicates the structure of the coefficients to use. \n\
  11           -c_str 0 => Setup for an isotropic material with constant coefficients. \n\
  12                          Parameters: \n\
  13                              -iso_E  : Youngs modulus \n\
  14                              -iso_nu : Poisson ratio \n\
  15           -c_str 1 => Setup for a step function in the material properties in x. \n\
  16                          Parameters: \n\
  17                               -step_E0  : Youngs modulus to the left of the step \n\
  18                               -step_nu0 : Poisson ratio to the left of the step \n\
  19                               -step_E1  : Youngs modulus to the right of the step \n\
  20                               -step_n1  : Poisson ratio to the right of the step \n\
  21                               -step_xc  : x coordinate of the step \n\
  22           -c_str 2 => Setup for a checkerboard material with alternating properties. \n\
  23                       Repeats the following pattern throughout the domain. For example with 4 materials specified, we would heve \n\
  24                       -------------------------\n\
  25                       |  D  |  A  |  B  |  C  |\n\
  26                       ------|-----|-----|------\n\
  27                       |  C  |  D  |  A  |  B  |\n\
  28                       ------|-----|-----|------\n\
  29                       |  B  |  C  |  D  |  A  |\n\
  30                       ------|-----|-----|------\n\
  31                       |  A  |  B  |  C  |  D  |\n\
  32                       -------------------------\n\
  33                       \n\
  34                          Parameters: \n\
  35                               -brick_E    : a comma seperated list of Young's modulii \n\
  36                               -brick_nu   : a comma seperated list of Poisson ratio's  \n\
  37                               -brick_span : the number of elements in x and y each brick will span \n\
  38           -c_str 3 => Setup for a sponge-like material with alternating properties. \n\
  39                       Repeats the following pattern throughout the domain \n\
  40                       -----------------------------\n\
  41                       |       [background]        |\n\
  42                       |          E0,nu0           |\n\
  43                       |     -----------------     |\n\
  44                       |     |  [inclusion]  |     |\n\
  45                       |     |    E1,nu1     |     |\n\
  46                       |     |               |     |\n\
  47                       |     | <---- w ----> |     |\n\
  48                       |     |               |     |\n\
  49                       |     |               |     |\n\
  50                       |     -----------------     |\n\
  51                       |                           |\n\
  52                       |                           |\n\
  53                       -----------------------------\n\
  54                       <--------  t + w + t ------->\n\
  55                       \n\
  56                          Parameters: \n\
  57                               -sponge_E0  : Youngs moduls of the surrounding material \n\
  58                               -sponge_E1  : Youngs moduls of the inclusio \n\
  59                               -sponge_nu0 : Poisson ratio of the surrounding material \n\
  60                               -sponge_nu1 : Poisson ratio of the inclusio \n\
  61                               -sponge_t   : the number of elements defining the border around each inclusion \n\
  62                               -sponge_w   : the number of elements in x and y each inclusion will span\n\
  63      -use_gp_coords : Evaluate the Youngs modulus, Poisson ratio and the body force at the global coordinates of the quadrature points.\n\
  64      By default, E, nu and the body force are evaulated at the element center and applied as a constant over the entire element.\n\
  65      -use_nonsymbc : Option to use non-symmetric boundary condition imposition. This choice will use less memory.";
  66
  67 /* Contributed by Dave May */
  68
  69 #include <petscksp.h>
  70 #include <petscdm.h>
  71 #include <petscdmda.h>
  72
  73 static PetscErrorCode DMDABCApplyCompression(DM,Mat,Vec);
  74 static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff);
  75
  76
  77 #define NSD            2 /* number of spatial dimensions */
  78 #define NODES_PER_EL   4 /* nodes per element */
  79 #define U_DOFS         2 /* degrees of freedom per displacement node */
  80 #define GAUSS_POINTS   4
  81
  82
  83
  84
  85 int KrylovMinimize(Mat A, Vec b, Vec x) {
  86
  87
  88   //Variables
  89
  90   PetscScalar  gamma, alpha, oldgamma, beta;
  91   PetscReal norm=20, Eprecision=1e-6, cgprec=1e-40;
  92   PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
  93   PetscErrorCode ierr;
  94   PetscScalar T1, T2;
  95   KSP ksp;
  96   PetscInt total=0;
  97   PetscInt size;
  98   PetscInt Istart,Iend;
  99   PetscInt i,its;
 100   Vec       x_old, residu;
 101   Mat S, AS;
 102   PetscScalar *array;
 103   PetscInt *ind_row;
 104   Vec Alpha, p, ss, vect, r, q, Ax;
 105
 106
 107   PetscInt first=1;
 108
 109   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
 110   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
 111   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 112
 113
 114
 115
 116   VecGetSize(b,&size);
 117
 118   ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);
 119
 120   PetscInt aa,bb;
 121   MatGetOwnershipRange(A,&aa,&bb);
 122
 123   //  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);
 124   //PetscSynchronizedFlush(PETSC_COMM_WORLD);
 125
 126
 127   ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
 128   ierr = MatSetSizes(S, bb-aa,  PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
 129   ierr = MatSetType(S, MATMPIDENSE); CHKERRQ(ierr);
 130   ierr = MatSetUp(S); CHKERRQ(ierr);
 131
 132   ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
 133
 134   ierr = VecCreate(PETSC_COMM_WORLD, &Alpha); CHKERRQ(ierr);
 135   ierr = VecSetSizes(Alpha, PETSC_DECIDE, ColS); CHKERRQ(ierr);
 136   ierr = VecSetFromOptions(Alpha); CHKERRQ(ierr);
 137   ierr = VecDuplicate(Alpha, &vect); CHKERRQ(ierr);
 138   ierr = VecDuplicate(Alpha, &p); CHKERRQ(ierr);
 139   ierr = VecDuplicate(Alpha, &ss); CHKERRQ(ierr);
 140   ierr = VecDuplicate(b, &r); CHKERRQ(ierr);
 141   ierr = VecDuplicate(b, &q); CHKERRQ(ierr);
 142   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
 143
 144   ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
 145   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
 146
 147
 148
 149   //indexes of row (these indexes are global)
 150   ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
 151   for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
 152
 153   //Initializations
 154   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
 155   ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
 156   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 157
 158
 159
 160   //GMRES WITH MINIMIZATION
 161   T1 = MPI_Wtime();
 162   ierr = KSPSetUp(ksp);CHKERRQ(ierr);
 163   while(giter<Emaxiter && norm>Eprecision ){
 164     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 165
 166
 167       //Solve
 168       ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
 169
 170       ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
 171       total += its; Iiter ++;
 172
 173
 174
 175       //Build S'
 176       ierr = VecGetArray(x, &array);
 177       ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
 178       VecRestoreArray(x, &array);
 179
 180
 181
 182       KSPGetResidualNorm(ksp,&norm);
 183
 184       /* //Error
 185       ierr = VecCopy(x, residu); CHKERRQ(ierr);
 186       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
 187       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
 188        */
 189
 190
 191       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
 192       ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 193
 194
 195     }
 196
 197
 198     //minimization step
 199     if( norm>Eprecision) {
 200
 201       MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
 202       MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
 203
 204
 205       //Build AS
 206       if(first) {
 207         MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
 208         first=0;
 209       }
 210       else
 211         MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
 212
 213
 214
 215
 216       //Minimization with CGLS
 217       MatMult(AS, Alpha, r);  VecAYPX(r, -1, b); //r_0 = b-AS*x_0
 218
 219
 220       MatMultTranspose(AS, r, p); //p_0 = AS'*r_0
 221
 222
 223
 224
 225       ierr = VecCopy(p, ss); CHKERRQ(ierr); //p_0 = s_0
 226       VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
 227       while(gamma>cgprec && iter<iterations){
 228         MatMult(AS, p, q); //q = AS*p
 229         VecNorm(q, NORM_2, &alpha); alpha = alpha * alpha; //alpha = norm2(q)^2
 230         alpha = gamma / alpha;
 231         VecAXPY(Alpha, alpha, p); //Alpha += alpha*p
 232         VecAXPY(r, -alpha, q); //r -= alpha*q
 233         MatMultTranspose(AS, r, ss); // ss = AS'*r
 234         oldgamma = gamma;
 235         VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
 236         beta = gamma / oldgamma;
 237         VecAYPX(p, beta, ss); //p = s+beta*p;
 238         iter ++;
 239       }
 240
 241       iter = 0; giter ++;
 242       //Minimizer
 243       MatMult(S, Alpha, x); //x = S*Alpha;
 244     }
 245
 246   }
 247   T2 = MPI_Wtime();
 248   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
 249   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
 250
 251   return 0;
 252
 253 }
 254
 255
 256
 257
 258
 259
 260
 261
 262
 263 int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 264
 265
 266   //Variables
 267
 268   PetscScalar  alpha, beta;
 269   PetscReal norm=20, Eprecision=1e-6, tol=1e-40;
 270   PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
 271   PetscErrorCode ierr;
 272   PetscScalar T1, T2;
 273   KSP ksp;
 274   PetscInt total=0;
 275   PetscInt size;
 276   PetscInt Istart,Iend;
 277   PetscInt i,its;
 278   Vec       x_old, residu;
 279   Mat S, AS;
 280   PetscScalar *array;
 281   PetscInt *ind_row;
 282   Vec  Ax;
 283   PetscScalar norm_ksp;
 284   Vec u,v,d,uu,vt,zero_long,zero_short,x_lsqr;
 285
 286   PetscInt first=1;
 287
 288   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
 289   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
 290   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 291
 292
 293
 294
 295   VecGetSize(b,&size);
 296
 297   ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);
 298
 299   PetscInt aa,bb;
 300   MatGetOwnershipRange(A,&aa,&bb);
 301
 302   //  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);
 303   //PetscSynchronizedFlush(PETSC_COMM_WORLD);
 304
 305
 306   ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
 307   ierr = MatSetSizes(S, bb-aa,  PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
 308   ierr = MatSetType(S, MATMPIDENSE); CHKERRQ(ierr);
 309   ierr = MatSetUp(S); CHKERRQ(ierr);
 310
 311   ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
 312
 313
 314
 315   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
 316
 317   ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
 318   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
 319
 320
 321   //long vector
 322   ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
 323
 324
 325   ierr = VecDuplicate(b,&uu);CHKERRQ(ierr);
 326   ierr = VecDuplicate(b,&zero_long);CHKERRQ(ierr);
 327   ierr = VecSet(zero_long,0);CHKERRQ(ierr);
 328
 329   //small vector
 330   ierr = VecCreate(PETSC_COMM_WORLD, &v); CHKERRQ(ierr);
 331   ierr = VecSetSizes(v, PETSC_DECIDE, ColS); CHKERRQ(ierr);
 332   ierr = VecSetFromOptions(v); CHKERRQ(ierr);
 333   ierr = VecDuplicate(v,&zero_short);CHKERRQ(ierr);
 334   ierr = VecSet(zero_short,0);CHKERRQ(ierr);
 335   ierr = VecDuplicate(v,&d);CHKERRQ(ierr);
 336   ierr = VecDuplicate(v,&vt);CHKERRQ(ierr);
 337   ierr = VecDuplicate(v,&x_lsqr);CHKERRQ(ierr);
 338
 339
 340   //indexes of row (these indexes are global)
 341   ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
 342   for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
 343
 344   //Initializations
 345   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
 346   ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
 347   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 348
 349
 350
 351
 352   //GMRES WITH MINIMIZATION
 353   T1 = MPI_Wtime();
 354   ierr = KSPSetUp(ksp);CHKERRQ(ierr);
 355   while(giter<Emaxiter && norm>Eprecision ){
 356     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 357
 358
 359       //Solve
 360       ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
 361
 362       ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
 363       total += its; Iiter ++;
 364
 365
 366
 367       //Build S'
 368       ierr = VecGetArray(x, &array);
 369       ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
 370       VecRestoreArray(x, &array);
 371
 372
 373
 374       KSPGetResidualNorm(ksp,&norm);
 375       /*
 376       //Error
 377       ierr = VecCopy(x, residu); CHKERRQ(ierr);
 378       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
 379       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
 380        */
 381
 382
 383       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
 384       ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 385
 386
 387     }
 388
 389
 390     //minimization step
 391     if( norm>Eprecision) {
 392
 393       MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
 394       MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
 395
 396
 397
 398
 399       //Build AS
 400       if(first) {
 401         MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
 402
 403         first=0;
 404       }
 405       else
 406         MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
 407
 408
 409
 410
 411
 412       //LSQR
 413       //LSQR
 414       //LSQR
 415
 416
 417
 418       PetscScalar n2b,tolb,normr,c,s,phibar,normar,norma,thet,rhot,rho,phi;
 419       PetscInt stag;
 420       tolb = tol * n2b;
 421       VecNorm(b, NORM_2, &n2b); //n2b = norm(b);
 422       ierr = VecCopy(b, u); CHKERRQ(ierr);   //u=b
 423       VecNorm(u, NORM_2, &beta); // beta=norm(u)
 424       normr=beta;
 425       if (beta != 0) {
 426         VecAYPX(u,1/beta,zero_long); //  u = u / beta;
 427       }
 428       c=1;
 429       s=0;
 430       phibar=beta;
 431       MatMultTranspose(AS, u, v); //v=A'*u
 432       ierr = VecSet(x_lsqr,0);CHKERRQ(ierr);
 433       VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
 434       if (alpha != 0) {
 435         VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
 436       }
 437       ierr = VecSet(d,0);CHKERRQ(ierr);
 438       normar = alpha * beta;
 439       norma=0;
 440       //stag=0;
 441       for(i=0;i<iterations;i++) {
 442         MatMult(AS, v, uu); //uu=A*v
 443         VecAYPX(u, -alpha, uu); //u = uu-alpha*u;
 444         VecNorm(u, NORM_2, &beta); // beta=norm(u)
 445         VecAYPX(u,1/beta,zero_long); //  u = u / beta;
 446         norma=sqrt(norma*norma+alpha*alpha+beta*beta); // norma = norm([norma alpha beta]);
 447         thet = - s * alpha;
 448         rhot = c * alpha;
 449         rho = sqrt(rhot*rhot + beta*beta);
 450         c = rhot / rho;
 451         s = - beta / rho;
 452         phi = c * phibar;
 453         if (phi == 0) {             // stagnation of the method
 454           stag = 1;
 455         }
 456         phibar = s * phibar;
 457         VecAYPX(d,-thet,v);       //d = (v - thet * d);
 458         VecAYPX(d,1/rho,zero_short);     //d=d/ rho;
 459
 460
 461
 462         VecAXPY(x_lsqr,phi,d);     // x_lsqr=x_lsqr+phi*d
 463         normr = abs(s) * normr;
 464         MatMultTranspose(AS, u, vt); //vt=A'*u;
 465         VecAYPX(v,-beta,vt);  // v = vt - beta * v;
 466         VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
 467         VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
 468         normar = alpha * abs( s * phi);
 469       }
 470
 471
 472
 473       iter = 0; giter ++;
 474       //Minimizer
 475       MatMult(S, x_lsqr, x); //x = S*x_small;
 476     }
 477
 478   }
 479   T2 = MPI_Wtime();
 480   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time LSQR : %g (s)\n", T2-T1); CHKERRQ(ierr);
 481   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations LSQR : %D\n", total); CHKERRQ(ierr);
 482
 483   return 0;
 484
 485 }
 486
 487
 488
 489
 490
 491
 492 /* cell based evaluation */
 493 typedef struct {
 494   PetscScalar E,nu,fx,fy;
 495 } Coefficients;
 496
 497 /* Gauss point based evaluation 8+4+4+4 = 20 */
 498 typedef struct {
 499   PetscScalar gp_coords[2*GAUSS_POINTS];
 500   PetscScalar E[GAUSS_POINTS];
 501   PetscScalar nu[GAUSS_POINTS];
 502   PetscScalar fx[GAUSS_POINTS];
 503   PetscScalar fy[GAUSS_POINTS];
 504 } GaussPointCoefficients;
 505
 506 typedef struct {
 507   PetscScalar ux_dof;
 508   PetscScalar uy_dof;
 509 } ElasticityDOF;
 510
 511
 512 /*
 513
 514  D = E/((1+nu)(1-2nu)) * [ 1-nu   nu        0     ]
 515                          [  nu   1-nu       0     ]
 516                          [  0     0   0.5*(1-2nu) ]
 517
 518  B = [ d_dx   0   ]
 519      [  0    d_dy ]
 520      [ d_dy  d_dx ]
 521
 522  */
 523
 524 /* FEM routines */
 525 /*
 526  Element: Local basis function ordering
 527  1-----2
 528  |     |
 529  |     |
 530  0-----3
 531  */
 532 static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
 533 {
 534   PetscScalar xi  = _xi[0];
 535   PetscScalar eta = _xi[1];
 536
 537   Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
 538   Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
 539   Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
 540   Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
 541 }
 542
 543 static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
 544 {
 545   PetscScalar xi  = _xi[0];
 546   PetscScalar eta = _xi[1];
 547
 548   GNi[0][0] = -0.25*(1.0-eta);
 549   GNi[0][1] = -0.25*(1.0+eta);
 550   GNi[0][2] =   0.25*(1.0+eta);
 551   GNi[0][3] =   0.25*(1.0-eta);
 552
 553   GNi[1][0] = -0.25*(1.0-xi);
 554   GNi[1][1] =   0.25*(1.0-xi);
 555   GNi[1][2] =   0.25*(1.0+xi);
 556   GNi[1][3] = -0.25*(1.0+xi);
 557 }
 558
 559 static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
 560 {
 561   PetscScalar J00,J01,J10,J11,J;
 562   PetscScalar iJ00,iJ01,iJ10,iJ11;
 563   PetscInt    i;
 564
 565   J00 = J01 = J10 = J11 = 0.0;
 566   for (i = 0; i < NODES_PER_EL; i++) {
 567     PetscScalar cx = coords[2*i+0];
 568     PetscScalar cy = coords[2*i+1];
 569
 570     J00 = J00+GNi[0][i]*cx;      /* J_xx = dx/dxi */
 571     J01 = J01+GNi[0][i]*cy;      /* J_xy = dy/dxi */
 572     J10 = J10+GNi[1][i]*cx;      /* J_yx = dx/deta */
 573     J11 = J11+GNi[1][i]*cy;      /* J_yy = dy/deta */
 574   }
 575   J = (J00*J11)-(J01*J10);
 576
 577   iJ00 =  J11/J;
 578   iJ01 = -J01/J;
 579   iJ10 = -J10/J;
 580   iJ11 =  J00/J;
 581
 582
 583   for (i = 0; i < NODES_PER_EL; i++) {
 584     GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
 585     GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
 586   }
 587
 588   if (det_J != NULL) *det_J = J;
 589 }
 590
 591 static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
 592 {
 593   *ngp         = 4;
 594   gp_xi[0][0]  = -0.57735026919;gp_xi[0][1] = -0.57735026919;
 595   gp_xi[1][0]  = -0.57735026919;gp_xi[1][1] =  0.57735026919;
 596   gp_xi[2][0]  =  0.57735026919;gp_xi[2][1] =  0.57735026919;
 597   gp_xi[3][0]  =  0.57735026919;gp_xi[3][1] = -0.57735026919;
 598   gp_weight[0] = 1.0;
 599   gp_weight[1] = 1.0;
 600   gp_weight[2] = 1.0;
 601   gp_weight[3] = 1.0;
 602 }
 603
 604
 605 /* procs to the left claim the ghost node as their element */
 606 #undef __FUNCT__
 607 #define __FUNCT__ "DMDAGetLocalElementSize"
 608 static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
 609 {
 610   PetscErrorCode ierr;
 611   PetscInt       m,n,p,M,N,P;
 612   PetscInt       sx,sy,sz;
 613
 614   PetscFunctionBeginUser;
 615   ierr = DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
 616   ierr = DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);CHKERRQ(ierr);
 617
 618   if (mxl != NULL) {
 619     *mxl = m;
 620     if ((sx+m) == M) *mxl = m-1;    /* last proc */
 621   }
 622   if (myl != NULL) {
 623     *myl = n;
 624     if ((sy+n) == N) *myl = n-1;  /* last proc */
 625   }
 626   if (mzl != NULL) {
 627     *mzl = p;
 628     if ((sz+p) == P) *mzl = p-1;  /* last proc */
 629   }
 630   PetscFunctionReturn(0);
 631 }
 632
 633 #undef __FUNCT__
 634 #define __FUNCT__ "DMDAGetElementCorners"
 635 static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
 636 {
 637   PetscErrorCode ierr;
 638   PetscInt       si,sj,sk;
 639
 640   PetscFunctionBeginUser;
 641   ierr = DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);CHKERRQ(ierr);
 642
 643   if (sx) {
 644     *sx = si;
 645     if (si != 0) *sx = si+1;
 646   }
 647   if (sy) {
 648     *sy = sj;
 649     if (sj != 0) *sy = sj+1;
 650   }
 651
 652   if (sk) {
 653     *sz = sk;
 654     if (sk != 0) *sz = sk+1;
 655   }
 656
 657   ierr = DMDAGetLocalElementSize(da,mx,my,mz);CHKERRQ(ierr);
 658   PetscFunctionReturn(0);
 659 }
 660
 661 #undef __FUNCT__
 662 #define __FUNCT__ "DMDAGetElementOwnershipRanges2d"
 663 static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
 664 {
 665   PetscErrorCode ierr;
 666   PetscMPIInt    rank;
 667   PetscInt       proc_I,proc_J;
 668   PetscInt       cpu_x,cpu_y;
 669   PetscInt       local_mx,local_my;
 670   Vec            vlx,vly;
 671   PetscInt       *LX,*LY,i;
 672   PetscScalar    *_a;
 673   Vec            V_SEQ;
 674   VecScatter     ctx;
 675
 676   PetscFunctionBeginUser;
 677   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 678
 679   DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
 680
 681   proc_J = rank/cpu_x;
 682   proc_I = rank-cpu_x*proc_J;
 683
 684   ierr = PetscMalloc(sizeof(PetscInt)*cpu_x,&LX);CHKERRQ(ierr);
 685   ierr = PetscMalloc(sizeof(PetscInt)*cpu_y,&LY);CHKERRQ(ierr);
 686
 687   ierr = DMDAGetLocalElementSize(da,&local_mx,&local_my,NULL);CHKERRQ(ierr);
 688   ierr = VecCreate(PETSC_COMM_WORLD,&vlx);CHKERRQ(ierr);
 689   ierr = VecSetSizes(vlx,PETSC_DECIDE,cpu_x);CHKERRQ(ierr);
 690   ierr = VecSetFromOptions(vlx);CHKERRQ(ierr);
 691
 692   ierr = VecCreate(PETSC_COMM_WORLD,&vly);CHKERRQ(ierr);
 693   ierr = VecSetSizes(vly,PETSC_DECIDE,cpu_y);CHKERRQ(ierr);
 694   ierr = VecSetFromOptions(vly);CHKERRQ(ierr);
 695
 696   ierr = VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);CHKERRQ(ierr);
 697   ierr = VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);CHKERRQ(ierr);
 698   ierr = VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);CHKERRQ(ierr);
 699   ierr = VecAssemblyBegin(vly);VecAssemblyEnd(vly);CHKERRQ(ierr);
 700
 701
 702
 703   ierr = VecScatterCreateToAll(vlx,&ctx,&V_SEQ);CHKERRQ(ierr);
 704   ierr = VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
 705   ierr = VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
 706   ierr = VecGetArray(V_SEQ,&_a);CHKERRQ(ierr);
 707   for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
 708   ierr = VecRestoreArray(V_SEQ,&_a);CHKERRQ(ierr);
 709   ierr = VecScatterDestroy(&ctx);CHKERRQ(ierr);
 710   ierr = VecDestroy(&V_SEQ);CHKERRQ(ierr);
 711
 712   ierr = VecScatterCreateToAll(vly,&ctx,&V_SEQ);CHKERRQ(ierr);
 713   ierr = VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
 714   ierr = VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
 715   ierr = VecGetArray(V_SEQ,&_a);CHKERRQ(ierr);
 716   for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
 717   ierr = VecRestoreArray(V_SEQ,&_a);CHKERRQ(ierr);
 718   ierr = VecScatterDestroy(&ctx);CHKERRQ(ierr);
 719   ierr = VecDestroy(&V_SEQ);CHKERRQ(ierr);
 720
 721   *_lx = LX;
 722   *_ly = LY;
 723
 724   ierr = VecDestroy(&vlx);CHKERRQ(ierr);
 725   ierr = VecDestroy(&vly);CHKERRQ(ierr);
 726   PetscFunctionReturn(0);
 727 }
 728
 729 #undef __FUNCT__
 730 #define __FUNCT__ "DMDACoordViewGnuplot2d"
 731 static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
 732 {
 733   DM             cda;
 734   Vec            coords;
 735   DMDACoor2d     **_coords;
 736   PetscInt       si,sj,nx,ny,i,j;
 737   FILE           *fp;
 738   char           fname[PETSC_MAX_PATH_LEN];
 739   PetscMPIInt    rank;
 740   PetscErrorCode ierr;
 741
 742   PetscFunctionBeginUser;
 743   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
 744   ierr = PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);CHKERRQ(ierr);
 745   ierr = PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);CHKERRQ(ierr);
 746   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
 747   ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);CHKERRQ(ierr);
 748
 749   ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
 750   ierr = DMGetCoordinatesLocal(da,&coords);CHKERRQ(ierr);
 751   ierr = DMDAVecGetArray(cda,coords,&_coords);CHKERRQ(ierr);
 752   ierr = DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);CHKERRQ(ierr);
 753   for (j = sj; j < sj+ny-1; j++) {
 754     for (i = si; i < si+nx-1; i++) {
 755       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));CHKERRQ(ierr);
 756       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i].x),PetscRealPart(_coords[j+1][i].y));CHKERRQ(ierr);
 757       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i+1].x),PetscRealPart(_coords[j+1][i+1].y));CHKERRQ(ierr);
 758       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i+1].x),PetscRealPart(_coords[j][i+1].y));CHKERRQ(ierr);
 759       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));CHKERRQ(ierr);
 760     }
 761   }
 762   ierr = DMDAVecRestoreArray(cda,coords,&_coords);CHKERRQ(ierr);
 763
 764   ierr = PetscFClose(PETSC_COMM_SELF,fp);CHKERRQ(ierr);
 765   PetscFunctionReturn(0);
 766 }
 767
 768 #undef __FUNCT__
 769 #define __FUNCT__ "DMDAViewGnuplot2d"
 770 static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
 771 {
 772   DM             cda;
 773   Vec            coords,local_fields;
 774   DMDACoor2d     **_coords;
 775   FILE           *fp;
 776   char           fname[PETSC_MAX_PATH_LEN];
 777   const char     *field_name;
 778   PetscMPIInt    rank;
 779   PetscInt       si,sj,nx,ny,i,j;
 780   PetscInt       n_dofs,d;
 781   PetscScalar    *_fields;
 782   PetscErrorCode ierr;
 783
 784   PetscFunctionBeginUser;
 785   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 786   ierr = PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);CHKERRQ(ierr);
 787   ierr = PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);CHKERRQ(ierr);
 788   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
 789
 790   ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);CHKERRQ(ierr);
 791   ierr = DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);CHKERRQ(ierr);
 792   ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");CHKERRQ(ierr);
 793   for (d = 0; d < n_dofs; d++) {
 794     ierr = DMDAGetFieldName(da,d,&field_name);CHKERRQ(ierr);
 795     ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);CHKERRQ(ierr);
 796   }
 797   ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");CHKERRQ(ierr);
 798
 799
 800   ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
 801   ierr = DMGetCoordinatesLocal(da,&coords);CHKERRQ(ierr);
 802   ierr = DMDAVecGetArray(cda,coords,&_coords);CHKERRQ(ierr);
 803   ierr = DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);CHKERRQ(ierr);
 804
 805   ierr = DMCreateLocalVector(da,&local_fields);CHKERRQ(ierr);
 806   ierr = DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);CHKERRQ(ierr);
 807   ierr = DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);CHKERRQ(ierr);
 808   ierr = VecGetArray(local_fields,&_fields);CHKERRQ(ierr);
 809
 810   for (j = sj; j < sj+ny; j++) {
 811     for (i = si; i < si+nx; i++) {
 812       PetscScalar coord_x,coord_y;
 813       PetscScalar field_d;
 814
 815       coord_x = _coords[j][i].x;
 816       coord_y = _coords[j][i].y;
 817
 818       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));CHKERRQ(ierr);
 819       for (d = 0; d < n_dofs; d++) {
 820         field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
 821         ierr    = PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",PetscRealPart(field_d));CHKERRQ(ierr);
 822       }
 823       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"\n");CHKERRQ(ierr);
 824     }
 825   }
 826   ierr = VecRestoreArray(local_fields,&_fields);CHKERRQ(ierr);
 827   ierr = VecDestroy(&local_fields);CHKERRQ(ierr);
 828
 829   ierr = DMDAVecRestoreArray(cda,coords,&_coords);CHKERRQ(ierr);
 830
 831   ierr = PetscFClose(PETSC_COMM_SELF,fp);CHKERRQ(ierr);
 832   PetscFunctionReturn(0);
 833 }
 834
 835 #undef __FUNCT__
 836 #define __FUNCT__ "DMDAViewCoefficientsGnuplot2d"
 837 static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
 838 {
 839   DM                     cda;
 840   Vec                    local_fields;
 841   FILE                   *fp;
 842   char                   fname[PETSC_MAX_PATH_LEN];
 843   const char             *field_name;
 844   PetscMPIInt            rank;
 845   PetscInt               si,sj,nx,ny,i,j,p;
 846   PetscInt               n_dofs,d;
 847   GaussPointCoefficients **_coefficients;
 848   PetscErrorCode         ierr;
 849
 850   PetscFunctionBeginUser;
 851   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
 852   ierr = PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);CHKERRQ(ierr);
 853   ierr = PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);CHKERRQ(ierr);
 854   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
 855
 856   ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);CHKERRQ(ierr);
 857   ierr = DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);CHKERRQ(ierr);
 858   ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");CHKERRQ(ierr);
 859   for (d = 0; d < n_dofs; d++) {
 860     ierr = DMDAGetFieldName(da,d,&field_name);CHKERRQ(ierr);
 861     ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);CHKERRQ(ierr);
 862   }
 863   ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");CHKERRQ(ierr);
 864
 865
 866   ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
 867   ierr = DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);CHKERRQ(ierr);
 868
 869   ierr = DMCreateLocalVector(da,&local_fields);CHKERRQ(ierr);
 870   ierr = DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);CHKERRQ(ierr);
 871   ierr = DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);CHKERRQ(ierr);
 872   ierr = DMDAVecGetArray(da,local_fields,&_coefficients);CHKERRQ(ierr);
 873
 874   for (j = sj; j < sj+ny; j++) {
 875     for (i = si; i < si+nx; i++) {
 876       PetscScalar coord_x,coord_y;
 877
 878       for (p = 0; p < GAUSS_POINTS; p++) {
 879         coord_x = _coefficients[j][i].gp_coords[2*p];
 880         coord_y = _coefficients[j][i].gp_coords[2*p+1];
 881
 882         ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));CHKERRQ(ierr);
 883
 884         ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e %1.6e",
 885                             PetscRealPart(_coefficients[j][i].E[p]),PetscRealPart(_coefficients[j][i].nu[p]),
 886                             PetscRealPart(_coefficients[j][i].fx[p]),PetscRealPart(_coefficients[j][i].fy[p]));CHKERRQ(ierr);
 887         ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"\n");CHKERRQ(ierr);
 888       }
 889     }
 890   }
 891   ierr = DMDAVecRestoreArray(da,local_fields,&_coefficients);CHKERRQ(ierr);
 892   ierr = VecDestroy(&local_fields);CHKERRQ(ierr);
 893
 894   ierr = PetscFClose(PETSC_COMM_SELF,fp);CHKERRQ(ierr);
 895   PetscFunctionReturn(0);
 896 }
 897
 898 static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar E[],PetscScalar nu[])
 899 {
 900   PetscInt    ngp;
 901   PetscScalar gp_xi[GAUSS_POINTS][2];
 902   PetscScalar gp_weight[GAUSS_POINTS];
 903   PetscInt    p,i,j,k,l;
 904   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
 905   PetscScalar J_p;
 906   PetscScalar B[3][U_DOFS*NODES_PER_EL];
 907   PetscScalar prop_E,prop_nu,factor,constit_D[3][3];
 908
 909   /* define quadrature rule */
 910   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
 911
 912   /* evaluate integral */
 913   for (p = 0; p < ngp; p++) {
 914     ConstructQ12D_GNi(gp_xi[p],GNi_p);
 915     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
 916
 917     for (i = 0; i < NODES_PER_EL; i++) {
 918       PetscScalar d_dx_i = GNx_p[0][i];
 919       PetscScalar d_dy_i = GNx_p[1][i];
 920
 921       B[0][2*i] = d_dx_i;  B[0][2*i+1] = 0.0;
 922       B[1][2*i] = 0.0;     B[1][2*i+1] = d_dy_i;
 923       B[2][2*i] = d_dy_i;  B[2][2*i+1] = d_dx_i;
 924     }
 925
 926     /* form D for the quadrature point */
 927     prop_E          = E[p];
 928     prop_nu         = nu[p];
 929     factor          = prop_E / ((1.0+prop_nu)*(1.0-2.0*prop_nu));
 930     constit_D[0][0] = 1.0-prop_nu;  constit_D[0][1] = prop_nu;      constit_D[0][2] = 0.0;
 931     constit_D[1][0] = prop_nu;      constit_D[1][1] = 1.0-prop_nu;  constit_D[1][2] = 0.0;
 932     constit_D[2][0] = 0.0;          constit_D[2][1] = 0.0;          constit_D[2][2] = 0.5*(1.0-2.0*prop_nu);
 933     for (i = 0; i < 3; i++) {
 934       for (j = 0; j < 3; j++) {
 935         constit_D[i][j] = factor * constit_D[i][j] * gp_weight[p] * J_p;
 936       }
 937     }
 938
 939     /* form Bt tildeD B */
 940     /*
 941      Ke_ij = Bt_ik . D_kl . B_lj
 942      = B_ki . D_kl . B_lj
 943      */
 944     for (i = 0; i < 8; i++) {
 945       for (j = 0; j < 8; j++) {
 946         for (k = 0; k < 3; k++) {
 947           for (l = 0; l < 3; l++) {
 948             Ke[8*i+j] = Ke[8*i+j] + B[k][i] * constit_D[k][l] * B[l][j];
 949           }
 950         }
 951       }
 952     }
 953
 954   } /* end quadrature */
 955 }
 956
 957 static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
 958 {
 959   PetscInt    ngp;
 960   PetscScalar gp_xi[GAUSS_POINTS][2];
 961   PetscScalar gp_weight[GAUSS_POINTS];
 962   PetscInt    p,i;
 963   PetscScalar Ni_p[NODES_PER_EL];
 964   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
 965   PetscScalar J_p,fac;
 966
 967   /* define quadrature rule */
 968   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
 969
 970   /* evaluate integral */
 971   for (p = 0; p < ngp; p++) {
 972     ConstructQ12D_Ni(gp_xi[p],Ni_p);
 973     ConstructQ12D_GNi(gp_xi[p],GNi_p);
 974     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
 975     fac = gp_weight[p]*J_p;
 976
 977     for (i = 0; i < NODES_PER_EL; i++) {
 978       Fe[NSD*i]   += fac*Ni_p[i]*fx[p];
 979       Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
 980     }
 981   }
 982 }
 983
 984 /*
 985  i,j are the element indices
 986  The unknown is a vector quantity.
 987  The s[].c is used to indicate the degree of freedom.
 988  */
 989 #undef __FUNCT__
 990 #define __FUNCT__ "DMDAGetElementEqnums_u"
 991 static PetscErrorCode DMDAGetElementEqnums_u(MatStencil s_u[],PetscInt i,PetscInt j)
 992 {
 993   PetscFunctionBeginUser;
 994   /* displacement */
 995   /* node 0 */
 996   s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0;          /* Ux0 */
 997   s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1;          /* Uy0 */
 998
 999   /* node 1 */
1000   s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0;        /* Ux1 */
1001   s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1;        /* Uy1 */
1002
1003   /* node 2 */
1004   s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0;      /* Ux2 */
1005   s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1;      /* Uy2 */
1006
1007   /* node 3 */
1008   s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0;        /* Ux3 */
1009   s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1;        /* Uy3 */
1010   PetscFunctionReturn(0);
1011 }
1012
1013 #undef __FUNCT__
1014 #define __FUNCT__ "GetElementCoords"
1015 static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
1016 {
1017   PetscFunctionBeginUser;
1018   /* get coords for the element */
1019   el_coords[NSD*0+0] = _coords[ej][ei].x;      el_coords[NSD*0+1] = _coords[ej][ei].y;
1020   el_coords[NSD*1+0] = _coords[ej+1][ei].x;    el_coords[NSD*1+1] = _coords[ej+1][ei].y;
1021   el_coords[NSD*2+0] = _coords[ej+1][ei+1].x;  el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
1022   el_coords[NSD*3+0] = _coords[ej][ei+1].x;    el_coords[NSD*3+1] = _coords[ej][ei+1].y;
1023   PetscFunctionReturn(0);
1024 }
1025
1026 #undef __FUNCT__
1027 #define __FUNCT__ "AssembleA_Elasticity"
1028 static PetscErrorCode AssembleA_Elasticity(Mat A,DM elas_da,DM properties_da,Vec properties)
1029 {
1030   DM                     cda;
1031   Vec                    coords;
1032   DMDACoor2d             **_coords;
1033   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
1034   PetscInt               sex,sey,mx,my;
1035   PetscInt               ei,ej;
1036   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
1037   PetscScalar            el_coords[NODES_PER_EL*NSD];
1038   Vec                    local_properties;
1039   GaussPointCoefficients **props;
1040   PetscScalar            *prop_E,*prop_nu;
1041   PetscErrorCode         ierr;
1042
1043   PetscFunctionBeginUser;
1044   /* setup for coords */
1045   ierr = DMGetCoordinateDM(elas_da,&cda);CHKERRQ(ierr);
1046   ierr = DMGetCoordinatesLocal(elas_da,&coords);CHKERRQ(ierr);
1047   ierr = DMDAVecGetArray(cda,coords,&_coords);CHKERRQ(ierr);
1048
1049   /* setup for coefficients */
1050   ierr = DMCreateLocalVector(properties_da,&local_properties);CHKERRQ(ierr);
1051   ierr = DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);CHKERRQ(ierr);
1052   ierr = DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);CHKERRQ(ierr);
1053   ierr = DMDAVecGetArray(properties_da,local_properties,&props);CHKERRQ(ierr);
1054
1055   ierr = DMDAGetElementCorners(elas_da,&sex,&sey,0,&mx,&my,0);CHKERRQ(ierr);
1056   for (ej = sey; ej < sey+my; ej++) {
1057     for (ei = sex; ei < sex+mx; ei++) {
1058       /* get coords for the element */
1059       GetElementCoords(_coords,ei,ej,el_coords);
1060
1061       /* get coefficients for the element */
1062       prop_E  = props[ej][ei].E;
1063       prop_nu = props[ej][ei].nu;
1064
1065       /* initialise element stiffness matrix */
1066       ierr = PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);CHKERRQ(ierr);
1067
1068       /* form element stiffness matrix */
1069       FormStressOperatorQ1(Ae,el_coords,prop_E,prop_nu);
1070
1071       /* insert element matrix into global matrix */
1072       ierr = DMDAGetElementEqnums_u(u_eqn,ei,ej);CHKERRQ(ierr);
1073       ierr = MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);CHKERRQ(ierr);
1074     }
1075   }
1076   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1077   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1078
1079   ierr = DMDAVecRestoreArray(cda,coords,&_coords);CHKERRQ(ierr);
1080
1081   ierr = DMDAVecRestoreArray(properties_da,local_properties,&props);CHKERRQ(ierr);
1082   ierr = VecDestroy(&local_properties);CHKERRQ(ierr);
1083   PetscFunctionReturn(0);
1084 }
1085
1086
1087 #undef __FUNCT__
1088 #define __FUNCT__ "DMDASetValuesLocalStencil_ADD_VALUES"
1089 static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F,MatStencil u_eqn[],PetscScalar Fe_u[])
1090 {
1091   PetscInt n;
1092
1093   PetscFunctionBeginUser;
1094   for (n = 0; n < 4; n++) {
1095     fields_F[u_eqn[2*n].j][u_eqn[2*n].i].ux_dof     = fields_F[u_eqn[2*n].j][u_eqn[2*n].i].ux_dof+Fe_u[2*n];
1096     fields_F[u_eqn[2*n+1].j][u_eqn[2*n+1].i].uy_dof = fields_F[u_eqn[2*n+1].j][u_eqn[2*n+1].i].uy_dof+Fe_u[2*n+1];
1097   }
1098   PetscFunctionReturn(0);
1099 }
1100
1101 #undef __FUNCT__
1102 #define __FUNCT__ "AssembleF_Elasticity"
1103 static PetscErrorCode AssembleF_Elasticity(Vec F,DM elas_da,DM properties_da,Vec properties)
1104 {
1105   DM                     cda;
1106   Vec                    coords;
1107   DMDACoor2d             **_coords;
1108   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
1109   PetscInt               sex,sey,mx,my;
1110   PetscInt               ei,ej;
1111   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
1112   PetscScalar            el_coords[NODES_PER_EL*NSD];
1113   Vec                    local_properties;
1114   GaussPointCoefficients **props;
1115   PetscScalar            *prop_fx,*prop_fy;
1116   Vec                    local_F;
1117   ElasticityDOF          **ff;
1118   PetscErrorCode         ierr;
1119
1120   PetscFunctionBeginUser;
1121   /* setup for coords */
1122   ierr = DMGetCoordinateDM(elas_da,&cda);CHKERRQ(ierr);
1123   ierr = DMGetCoordinatesLocal(elas_da,&coords);CHKERRQ(ierr);
1124   ierr = DMDAVecGetArray(cda,coords,&_coords);CHKERRQ(ierr);
1125
1126   /* setup for coefficients */
1127   ierr = DMGetLocalVector(properties_da,&local_properties);CHKERRQ(ierr);
1128   ierr = DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);CHKERRQ(ierr);
1129   ierr = DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);CHKERRQ(ierr);
1130   ierr = DMDAVecGetArray(properties_da,local_properties,&props);CHKERRQ(ierr);
1131
1132   /* get acces to the vector */
1133   ierr = DMGetLocalVector(elas_da,&local_F);CHKERRQ(ierr);
1134   ierr = VecZeroEntries(local_F);CHKERRQ(ierr);
1135   ierr = DMDAVecGetArray(elas_da,local_F,&ff);CHKERRQ(ierr);
1136
1137
1138   ierr = DMDAGetElementCorners(elas_da,&sex,&sey,0,&mx,&my,0);CHKERRQ(ierr);
1139   for (ej = sey; ej < sey+my; ej++) {
1140     for (ei = sex; ei < sex+mx; ei++) {
1141       /* get coords for the element */
1142       GetElementCoords(_coords,ei,ej,el_coords);
1143
1144       /* get coefficients for the element */
1145       prop_fx = props[ej][ei].fx;
1146       prop_fy = props[ej][ei].fy;
1147
1148       /* initialise element stiffness matrix */
1149       ierr = PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);CHKERRQ(ierr);
1150
1151       /* form element stiffness matrix */
1152       FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
1153
1154       /* insert element matrix into global matrix */
1155       ierr = DMDAGetElementEqnums_u(u_eqn,ei,ej);CHKERRQ(ierr);
1156
1157       ierr = DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,Fe);CHKERRQ(ierr);
1158     }
1159   }
1160
1161   ierr = DMDAVecRestoreArray(elas_da,local_F,&ff);CHKERRQ(ierr);
1162   ierr = DMLocalToGlobalBegin(elas_da,local_F,ADD_VALUES,F);CHKERRQ(ierr);
1163   ierr = DMLocalToGlobalEnd(elas_da,local_F,ADD_VALUES,F);CHKERRQ(ierr);
1164   ierr = DMRestoreLocalVector(elas_da,&local_F);CHKERRQ(ierr);
1165
1166   ierr = DMDAVecRestoreArray(cda,coords,&_coords);CHKERRQ(ierr);
1167
1168   ierr = DMDAVecRestoreArray(properties_da,local_properties,&props);CHKERRQ(ierr);
1169   ierr = DMRestoreLocalVector(properties_da,&local_properties);CHKERRQ(ierr);
1170   PetscFunctionReturn(0);
1171 }
1172
1173 #undef __FUNCT__
1174 #define __FUNCT__ "solve_elasticity_2d"
1175 static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
1176 {
1177   DM                     elas_da,da_prop;
1178   PetscInt               u_dof,dof,stencil_width;
1179   Mat                    A;
1180   PetscInt               mxl,myl;
1181   DM                     prop_cda,vel_cda;
1182   Vec                    prop_coords,vel_coords;
1183   PetscInt               si,sj,nx,ny,i,j,p;
1184   Vec                    f,X;
1185   PetscInt               prop_dof,prop_stencil_width;
1186   Vec                    properties,l_properties;
1187   MatNullSpace           matnull;
1188   PetscReal              dx,dy;
1189   PetscInt               M,N;
1190   DMDACoor2d             **_prop_coords,**_vel_coords;
1191   GaussPointCoefficients **element_props;
1192   KSP                    ksp_E;
1193   PetscInt               coefficient_structure = 0;
1194   PetscInt               cpu_x,cpu_y,*lx = NULL,*ly = NULL;
1195   PetscBool              use_gp_coords = PETSC_FALSE;
1196   PetscBool              use_nonsymbc  = PETSC_FALSE;
1197   PetscBool              no_view       = PETSC_FALSE;
1198   PetscBool              flg;
1199   PetscErrorCode         ierr;
1200
1201   PetscFunctionBeginUser;
1202   /* Generate the da for velocity and pressure */
1203   /*
1204    We use Q1 elements for the temperature.
1205    FEM has a 9-point stencil (BOX) or connectivity pattern
1206    Num nodes in each direction is mx+1, my+1
1207    */
1208   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1209   dof           = u_dof;
1210   stencil_width = 1;
1211   ierr          = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1212                                mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&elas_da);CHKERRQ(ierr);
1213   ierr = DMDASetFieldName(elas_da,0,"Ux");CHKERRQ(ierr);
1214   ierr = DMDASetFieldName(elas_da,1,"Uy");CHKERRQ(ierr);
1215
1216   /* unit box [0,1] x [0,1] */
1217   ierr = DMDASetUniformCoordinates(elas_da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);
1218
1219
1220   /* Generate element properties, we will assume all material properties are constant over the element */
1221   /* local number of elements */
1222   ierr = DMDAGetLocalElementSize(elas_da,&mxl,&myl,NULL);CHKERRQ(ierr);
1223
1224   /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1225   ierr = DMDAGetInfo(elas_da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);CHKERRQ(ierr);
1226   ierr = DMDAGetElementOwnershipRanges2d(elas_da,&lx,&ly);CHKERRQ(ierr);
1227
1228   prop_dof           = (PetscInt)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1229   prop_stencil_width = 0;
1230   ierr               = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1231                                     mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);CHKERRQ(ierr);
1232   ierr = PetscFree(lx);CHKERRQ(ierr);
1233   ierr = PetscFree(ly);CHKERRQ(ierr);
1234
1235   /* define centroid positions */
1236   ierr = DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
1237   dx   = 1.0/((PetscReal)(M));
1238   dy   = 1.0/((PetscReal)(N));
1239
1240   ierr = DMDASetUniformCoordinates(da_prop,0.0+0.5*dx,1.0-0.5*dx,0.0+0.5*dy,1.0-0.5*dy,0.0,1.0);CHKERRQ(ierr);
1241
1242   /* define coefficients */
1243   ierr = PetscOptionsGetInt(NULL,"-c_str",&coefficient_structure,NULL);CHKERRQ(ierr);
1244
1245   ierr = DMCreateGlobalVector(da_prop,&properties);CHKERRQ(ierr);
1246   ierr = DMCreateLocalVector(da_prop,&l_properties);CHKERRQ(ierr);
1247   ierr = DMDAVecGetArray(da_prop,l_properties,&element_props);CHKERRQ(ierr);
1248
1249   ierr = DMGetCoordinateDM(da_prop,&prop_cda);CHKERRQ(ierr);
1250   ierr = DMGetCoordinatesLocal(da_prop,&prop_coords);CHKERRQ(ierr);
1251   ierr = DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);CHKERRQ(ierr);
1252
1253   ierr = DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);CHKERRQ(ierr);
1254
1255   ierr = DMGetCoordinateDM(elas_da,&vel_cda);CHKERRQ(ierr);
1256   ierr = DMGetCoordinatesLocal(elas_da,&vel_coords);CHKERRQ(ierr);
1257   ierr = DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);CHKERRQ(ierr);
1258
1259
1260   /* interpolate the coordinates */
1261   for (j = sj; j < sj+ny; j++) {
1262     for (i = si; i < si+nx; i++) {
1263       PetscInt    ngp;
1264       PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1265       PetscScalar el_coords[8];
1266
1267       ierr = GetElementCoords(_vel_coords,i,j,el_coords);CHKERRQ(ierr);
1268       ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1269
1270       for (p = 0; p < GAUSS_POINTS; p++) {
1271         PetscScalar gp_x,gp_y;
1272         PetscInt    n;
1273         PetscScalar xi_p[2],Ni_p[4];
1274
1275         xi_p[0] = gp_xi[p][0];
1276         xi_p[1] = gp_xi[p][1];
1277         ConstructQ12D_Ni(xi_p,Ni_p);
1278
1279         gp_x = 0.0;
1280         gp_y = 0.0;
1281         for (n = 0; n < NODES_PER_EL; n++) {
1282           gp_x = gp_x+Ni_p[n]*el_coords[2*n];
1283           gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1284         }
1285         element_props[j][i].gp_coords[2*p]   = gp_x;
1286         element_props[j][i].gp_coords[2*p+1] = gp_y;
1287       }
1288     }
1289   }
1290
1291   /* define the coefficients */
1292   ierr = PetscOptionsGetBool(NULL,"-use_gp_coords",&use_gp_coords,&flg);CHKERRQ(ierr);
1293
1294   for (j = sj; j < sj+ny; j++) {
1295     for (i = si; i < si+nx; i++) {
1296       PetscScalar              centroid_x = _prop_coords[j][i].x; /* centroids of cell */
1297       PetscScalar              centroid_y = _prop_coords[j][i].y;
1298       PETSC_UNUSED PetscScalar coord_x,coord_y;
1299
1300
1301       if (coefficient_structure == 0) { /* isotropic */
1302         PetscScalar opts_E,opts_nu;
1303
1304         opts_E  = 1.0;
1305         opts_nu = 0.33;
1306         ierr    = PetscOptionsGetScalar(NULL,"-iso_E",&opts_E,&flg);CHKERRQ(ierr);
1307         ierr    = PetscOptionsGetScalar(NULL,"-iso_nu",&opts_nu,&flg);CHKERRQ(ierr);
1308
1309         for (p = 0; p < GAUSS_POINTS; p++) {
1310           element_props[j][i].E[p]  = opts_E;
1311           element_props[j][i].nu[p] = opts_nu;
1312
1313           element_props[j][i].fx[p] = 0.0;
1314           element_props[j][i].fy[p] = 0.0;
1315         }
1316       } else if (coefficient_structure == 1) { /* step */
1317         PetscScalar opts_E0,opts_nu0,opts_xc;
1318         PetscScalar opts_E1,opts_nu1;
1319
1320         opts_E0  = opts_E1  = 1.0;
1321         opts_nu0 = opts_nu1 = 0.333;
1322         opts_xc  = 0.5;
1323         ierr     = PetscOptionsGetScalar(NULL,"-step_E0",&opts_E0,&flg);CHKERRQ(ierr);
1324         ierr     = PetscOptionsGetScalar(NULL,"-step_nu0",&opts_nu0,&flg);CHKERRQ(ierr);
1325         ierr     = PetscOptionsGetScalar(NULL,"-step_E1",&opts_E1,&flg);CHKERRQ(ierr);
1326         ierr     = PetscOptionsGetScalar(NULL,"-step_nu1",&opts_nu1,&flg);CHKERRQ(ierr);
1327         ierr     = PetscOptionsGetScalar(NULL,"-step_xc",&opts_xc,&flg);CHKERRQ(ierr);
1328
1329         for (p = 0; p < GAUSS_POINTS; p++) {
1330           coord_x = centroid_x;
1331           coord_y = centroid_y;
1332           if (use_gp_coords) {
1333             coord_x = element_props[j][i].gp_coords[2*p];
1334             coord_y = element_props[j][i].gp_coords[2*p+1];
1335           }
1336
1337           element_props[j][i].E[p]  = opts_E0;
1338           element_props[j][i].nu[p] = opts_nu0;
1339           if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
1340             element_props[j][i].E[p]  = opts_E1;
1341             element_props[j][i].nu[p] = opts_nu1;
1342           }
1343
1344           element_props[j][i].fx[p] = 0.0;
1345           element_props[j][i].fy[p] = 0.0;
1346         }
1347       } else if (coefficient_structure == 2) { /* brick */
1348         PetscReal values_E[10];
1349         PetscReal values_nu[10];
1350         PetscInt  nbricks,maxnbricks;
1351         PetscInt  index,span;
1352         PetscInt  jj;
1353
1354         flg        = PETSC_FALSE;
1355         maxnbricks = 10;
1356         ierr       = PetscOptionsGetRealArray(NULL, "-brick_E",values_E,&maxnbricks,&flg);CHKERRQ(ierr);
1357         nbricks    = maxnbricks;
1358         if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of E values for each brick");CHKERRQ(ierr);
1359
1360         flg        = PETSC_FALSE;
1361         maxnbricks = 10;
1362         ierr       = PetscOptionsGetRealArray(NULL, "-brick_nu",values_nu,&maxnbricks,&flg);CHKERRQ(ierr);
1363         if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of nu values for each brick");CHKERRQ(ierr);
1364         if (maxnbricks != nbricks) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply equal numbers of values for E and nu");CHKERRQ(ierr);
1365
1366         span = 1;
1367         ierr = PetscOptionsGetInt(NULL,"-brick_span",&span,&flg);CHKERRQ(ierr);
1368
1369         /* cycle through the indices so that no two material properties are repeated in lines of x or y */
1370         jj    = (j/span)%nbricks;
1371         index = (jj+i/span)%nbricks;
1372         /*printf("j=%d: index = %d \n", j,index); */
1373
1374         for (p = 0; p < GAUSS_POINTS; p++) {
1375           element_props[j][i].E[p]  = values_E[index];
1376           element_props[j][i].nu[p] = values_nu[index];
1377         }
1378       } else if (coefficient_structure == 3) { /* sponge */
1379         PetscScalar opts_E0,opts_nu0;
1380         PetscScalar opts_E1,opts_nu1;
1381         PetscInt    opts_t,opts_w;
1382         PetscInt    ii,jj,ci,cj;
1383
1384         opts_E0  = opts_E1  = 1.0;
1385         opts_nu0 = opts_nu1 = 0.333;
1386         ierr     = PetscOptionsGetScalar(NULL,"-sponge_E0",&opts_E0,&flg);CHKERRQ(ierr);
1387         ierr     = PetscOptionsGetScalar(NULL,"-sponge_nu0",&opts_nu0,&flg);CHKERRQ(ierr);
1388         ierr     = PetscOptionsGetScalar(NULL,"-sponge_E1",&opts_E1,&flg);CHKERRQ(ierr);
1389         ierr     = PetscOptionsGetScalar(NULL,"-sponge_nu1",&opts_nu1,&flg);CHKERRQ(ierr);
1390
1391         opts_t = opts_w = 1;
1392         ierr   = PetscOptionsGetInt(NULL,"-sponge_t",&opts_t,&flg);CHKERRQ(ierr);
1393         ierr   = PetscOptionsGetInt(NULL,"-sponge_w",&opts_w,&flg);CHKERRQ(ierr);
1394
1395         ii = (i)/(opts_t+opts_w+opts_t);
1396         jj = (j)/(opts_t+opts_w+opts_t);
1397
1398         ci = i - ii*(opts_t+opts_w+opts_t);
1399         cj = j - jj*(opts_t+opts_w+opts_t);
1400
1401         for (p = 0; p < GAUSS_POINTS; p++) {
1402           element_props[j][i].E[p]  = opts_E0;
1403           element_props[j][i].nu[p] = opts_nu0;
1404         }
1405         if ((ci >= opts_t) && (ci < opts_t+opts_w)) {
1406           if ((cj >= opts_t) && (cj < opts_t+opts_w)) {
1407             for (p = 0; p < GAUSS_POINTS; p++) {
1408               element_props[j][i].E[p]  = opts_E1;
1409               element_props[j][i].nu[p] = opts_nu1;
1410             }
1411           }
1412         }
1413
1414       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1415     }
1416   }
1417   ierr = DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);CHKERRQ(ierr);
1418
1419   ierr = DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);CHKERRQ(ierr);
1420
1421   ierr = DMDAVecRestoreArray(da_prop,l_properties,&element_props);CHKERRQ(ierr);
1422   ierr = DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);CHKERRQ(ierr);
1423   ierr = DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);CHKERRQ(ierr);
1424
1425   /*  ierr = PetscOptionsGetBool(NULL,"-no_view",&no_view,NULL);CHKERRQ(ierr);
1426   if (!no_view) {
1427     ierr = DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for elasticity eqn.","properties");CHKERRQ(ierr);
1428     ierr = DMDACoordViewGnuplot2d(elas_da,"mesh");CHKERRQ(ierr);
1429    }*/
1430
1431   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1432   ierr = DMSetMatType(elas_da,MATAIJ);CHKERRQ(ierr);
1433   ierr = DMCreateMatrix(elas_da,&A);CHKERRQ(ierr);
1434   ierr = DMGetCoordinates(elas_da,&vel_coords);CHKERRQ(ierr);
1435   ierr = MatNullSpaceCreateRigidBody(vel_coords,&matnull);CHKERRQ(ierr);
1436   ierr = MatSetNearNullSpace(A,matnull);CHKERRQ(ierr);
1437   ierr = MatNullSpaceDestroy(&matnull);CHKERRQ(ierr);
1438   ierr = MatGetVecs(A,&f,&X);CHKERRQ(ierr);
1439
1440   /* assemble A11 */
1441   ierr = MatZeroEntries(A);CHKERRQ(ierr);
1442   ierr = VecZeroEntries(f);CHKERRQ(ierr);
1443
1444   ierr = AssembleA_Elasticity(A,elas_da,da_prop,properties);CHKERRQ(ierr);
1445   /* build force vector */
1446   ierr = AssembleF_Elasticity(f,elas_da,da_prop,properties);CHKERRQ(ierr);
1447
1448
1449   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_E);CHKERRQ(ierr);
1450   // ierr = KSPSetOptionsPrefix(ksp_E,"elas_");CHKERRQ(ierr);  /* elasticity */
1451
1452   //ierr = PetscOptionsGetBool(NULL,"-use_nonsymbc",&use_nonsymbc,&flg);CHKERRQ(ierr);
1453   /* solve */
1454   if (!use_nonsymbc) {
1455     Mat        AA;
1456     Vec        ff,XX;
1457     IS         is;
1458     VecScatter scat;
1459
1460     ierr = DMDABCApplySymmetricCompression(elas_da,A,f,&is,&AA,&ff);CHKERRQ(ierr);
1461     ierr = VecDuplicate(ff,&XX);CHKERRQ(ierr);
1462
1463     ierr = KSPSetOperators(ksp_E,AA,AA);CHKERRQ(ierr);
1464     ierr = KSPSetFromOptions(ksp_E);CHKERRQ(ierr);
1465
1466
1467    PetscScalar T1,T2;
1468     ierr = KSPSetTolerances(ksp_E, 1e-7, 1e-7, PETSC_DEFAULT, 50000000); CHKERRQ(ierr);
1469
1470
1471     PC             pc;
1472     KSPGetPC(ksp_E, &pc);
1473     PCType         type;
1474     PCGetType(pc, &type);
1475
1476     PetscPrintf(PETSC_COMM_WORLD, "PC TYPE %s  \n", type);
1477     KSPGetType(ksp_E,&type);
1478     PetscPrintf(PETSC_COMM_WORLD, "SOLVER TYPE %s  \n", type);
1479
1480
1481     T1 = MPI_Wtime();
1482  ierr = KSPSetUp(ksp_E);CHKERRQ(ierr);
1483     ierr = KSPSolve(ksp_E,ff,XX);CHKERRQ(ierr);
1484     T2 = MPI_Wtime();
1485
1486     Mat A;
1487     Vec sol;
1488     PetscScalar norm;
1489     Vec X2;
1490
1491     VecDuplicate(ff,&sol);
1492     VecDuplicate(ff,&X2);
1493     VecCopy(ff,X2);
1494
1495
1496     KSPGetOperators(ksp_E,&A,NULL);
1497     MatMult(A,XX,sol);
1498     VecAXPY(sol,-1,ff);
1499     VecNorm(sol, NORM_2, &norm);
1500
1501
1502     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Norm of error : %g\n", (double)norm); CHKERRQ(ierr);
1503     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
1504     PetscInt total;
1505     ierr = KSPGetIterationNumber(ksp_E, &total); CHKERRQ(ierr);
1506     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
1507
1508
1509
1510
1511     KrylovMinimize(A, ff, X2);
1512     MatMult(A,X2,sol);
1513
1514
1515     VecAXPY(sol,-1,ff);
1516     VecNorm(sol, NORM_2, &norm);
1517     ierr = PetscPrintf(PETSC_COMM_WORLD, "Error2 %g\n",norm);
1518
1519
1520
1521
1522     VecCopy(ff,X2);
1523     KrylovMinimizeLSQR(A, ff, X2);
1524     MatMult(A,X2,sol);
1525
1526
1527     VecAXPY(sol,-1,ff);
1528     VecNorm(sol, NORM_2, &norm);
1529     ierr = PetscPrintf(PETSC_COMM_WORLD, "ErrorLSQR %g\n",norm);
1530
1531
1532
1533     /* push XX back into X */
1534     ierr = DMDABCApplyCompression(elas_da,NULL,X);CHKERRQ(ierr);
1535
1536     ierr = VecScatterCreate(XX,NULL,X,is,&scat);CHKERRQ(ierr);
1537     ierr = VecScatterBegin(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1538     ierr = VecScatterEnd(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1539     ierr = VecScatterDestroy(&scat);CHKERRQ(ierr);
1540
1541     ierr = MatDestroy(&AA);CHKERRQ(ierr);
1542     ierr = VecDestroy(&ff);CHKERRQ(ierr);
1543     ierr = VecDestroy(&XX);CHKERRQ(ierr);
1544     ierr = ISDestroy(&is);CHKERRQ(ierr);
1545   } else {
1546     ierr = DMDABCApplyCompression(elas_da,A,f);CHKERRQ(ierr);
1547
1548     ierr = KSPSetOperators(ksp_E,A,A);CHKERRQ(ierr);
1549     ierr = KSPSetFromOptions(ksp_E);CHKERRQ(ierr);
1550
1551     ierr = KSPSolve(ksp_E,f,X);CHKERRQ(ierr);
1552   }
1553
1554   //  if (!no_view) {ierr = DMDAViewGnuplot2d(elas_da,X,"Displacement solution for elasticity eqn.","X");CHKERRQ(ierr);}
1555   ierr = KSPDestroy(&ksp_E);CHKERRQ(ierr);
1556
1557   ierr = VecDestroy(&X);CHKERRQ(ierr);
1558   ierr = VecDestroy(&f);CHKERRQ(ierr);
1559   ierr = MatDestroy(&A);CHKERRQ(ierr);
1560
1561   ierr = DMDestroy(&elas_da);CHKERRQ(ierr);
1562   ierr = DMDestroy(&da_prop);CHKERRQ(ierr);
1563
1564   ierr = VecDestroy(&properties);CHKERRQ(ierr);
1565   ierr = VecDestroy(&l_properties);CHKERRQ(ierr);
1566   PetscFunctionReturn(0);
1567 }
1568
1569 #undef __FUNCT__
1570 #define __FUNCT__ "main"
1571 int main(int argc,char **args)
1572 {
1573   PetscErrorCode ierr;
1574   PetscInt       mx,my;
1575
1576   ierr = PetscInitialize(&argc,&args,(char*)0,help);CHKERRQ(ierr);
1577
1578   mx   = my = 10;
1579   ierr = PetscOptionsGetInt(NULL,"-mx",&mx,NULL);CHKERRQ(ierr);
1580   ierr = PetscOptionsGetInt(NULL,"-my",&my,NULL);CHKERRQ(ierr);
1581
1582
1583   PetscMPIInt size;
1584   MPI_Comm_size(PETSC_COMM_WORLD,&size);
1585   PetscPrintf(PETSC_COMM_WORLD,"Number of processors = %d\n",size);
1586
1587
1588   ierr = solve_elasticity_2d(mx,my);CHKERRQ(ierr);
1589
1590   ierr = PetscFinalize();CHKERRQ(ierr);
1591   return 0;
1592 }
1593
1594 /* -------------------------- helpers for boundary conditions -------------------------------- */
1595
1596 #undef __FUNCT__
1597 #define __FUNCT__ "BCApply_EAST"
1598 static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1599 {
1600   DM                     cda;
1601   Vec                    coords;
1602   PetscInt               si,sj,nx,ny,i,j;
1603   PetscInt               M,N;
1604   DMDACoor2d             **_coords;
1605   const PetscInt         *g_idx;
1606   PetscInt               *bc_global_ids;
1607   PetscScalar            *bc_vals;
1608   PetscInt               nbcs;
1609   PetscInt               n_dofs;
1610   PetscErrorCode         ierr;
1611   ISLocalToGlobalMapping ltogm;
1612
1613   PetscFunctionBeginUser;
1614   /* enforce bc's */
1615   ierr = DMGetLocalToGlobalMapping(da,&ltogm);CHKERRQ(ierr);
1616   ierr = ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);CHKERRQ(ierr);
1617
1618   ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
1619   ierr = DMGetCoordinatesLocal(da,&coords);CHKERRQ(ierr);
1620   ierr = DMDAVecGetArray(cda,coords,&_coords);CHKERRQ(ierr);
1621   ierr = DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);CHKERRQ(ierr);
1622   ierr = DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);CHKERRQ(ierr);
1623
1624   /* --- */
1625
1626   ierr = PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);CHKERRQ(ierr);
1627   ierr = PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);CHKERRQ(ierr);
1628
1629   /* init the entries to -1 so VecSetValues will ignore them */
1630   for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1631
1632   i = nx-1;
1633   for (j = 0; j < ny; j++) {
1634     PetscInt                 local_id;
1635     PETSC_UNUSED PetscScalar coordx,coordy;
1636
1637     local_id = i+j*nx;
1638
1639     bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1640
1641     coordx = _coords[j+sj][i+si].x;
1642     coordy = _coords[j+sj][i+si].y;
1643
1644     bc_vals[j] =  bc_val;
1645   }
1646   ierr = ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);CHKERRQ(ierr);
1647   nbcs = 0;
1648   if ((si+nx) == (M)) nbcs = ny;
1649
1650   if (b) {
1651     ierr = VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);CHKERRQ(ierr);
1652     ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
1653     ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
1654   }
1655   if (A) {
1656     ierr = MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);CHKERRQ(ierr);
1657   }
1658
1659   ierr = PetscFree(bc_vals);CHKERRQ(ierr);
1660   ierr = PetscFree(bc_global_ids);CHKERRQ(ierr);
1661
1662   ierr = DMDAVecRestoreArray(cda,coords,&_coords);CHKERRQ(ierr);
1663   PetscFunctionReturn(0);
1664 }
1665
1666 #undef __FUNCT__
1667 #define __FUNCT__ "BCApply_WEST"
1668 static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1669 {
1670   DM                     cda;
1671   Vec                    coords;
1672   PetscInt               si,sj,nx,ny,i,j;
1673   PetscInt               M,N;
1674   DMDACoor2d             **_coords;
1675   const PetscInt         *g_idx;
1676   PetscInt               *bc_global_ids;
1677   PetscScalar            *bc_vals;
1678   PetscInt               nbcs;
1679   PetscInt               n_dofs;
1680   PetscErrorCode         ierr;
1681   ISLocalToGlobalMapping ltogm;
1682
1683   PetscFunctionBeginUser;
1684   /* enforce bc's */
1685   ierr = DMGetLocalToGlobalMapping(da,&ltogm);CHKERRQ(ierr);
1686   ierr = ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);CHKERRQ(ierr);
1687
1688   ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
1689   ierr = DMGetCoordinatesLocal(da,&coords);CHKERRQ(ierr);
1690   ierr = DMDAVecGetArray(cda,coords,&_coords);CHKERRQ(ierr);
1691   ierr = DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);CHKERRQ(ierr);
1692   ierr = DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);CHKERRQ(ierr);
1693
1694   /* --- */
1695
1696   ierr = PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);CHKERRQ(ierr);
1697   ierr = PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);CHKERRQ(ierr);
1698
1699   /* init the entries to -1 so VecSetValues will ignore them */
1700   for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1701
1702   i = 0;
1703   for (j = 0; j < ny; j++) {
1704     PetscInt                 local_id;
1705     PETSC_UNUSED PetscScalar coordx,coordy;
1706
1707     local_id = i+j*nx;
1708
1709     bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1710
1711     coordx = _coords[j+sj][i+si].x;
1712     coordy = _coords[j+sj][i+si].y;
1713
1714     bc_vals[j] =  bc_val;
1715   }
1716   ierr = ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);CHKERRQ(ierr);
1717   nbcs = 0;
1718   if (si == 0) nbcs = ny;
1719
1720   if (b) {
1721     ierr = VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);CHKERRQ(ierr);
1722     ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
1723     ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
1724   }
1725   if (A) {
1726     ierr = MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);CHKERRQ(ierr);
1727   }
1728
1729   ierr = PetscFree(bc_vals);CHKERRQ(ierr);
1730   ierr = PetscFree(bc_global_ids);CHKERRQ(ierr);
1731
1732   ierr = DMDAVecRestoreArray(cda,coords,&_coords);CHKERRQ(ierr);
1733   PetscFunctionReturn(0);
1734 }
1735
1736 #undef __FUNCT__
1737 #define __FUNCT__ "DMDABCApplyCompression"
1738 static PetscErrorCode DMDABCApplyCompression(DM elas_da,Mat A,Vec f)
1739 {
1740   PetscErrorCode ierr;
1741
1742   PetscFunctionBeginUser;
1743   ierr = BCApply_EAST(elas_da,0,-1.0,A,f);CHKERRQ(ierr);
1744   ierr = BCApply_EAST(elas_da,1, 0.0,A,f);CHKERRQ(ierr);
1745   ierr = BCApply_WEST(elas_da,0,1.0,A,f);CHKERRQ(ierr);
1746   ierr = BCApply_WEST(elas_da,1,0.0,A,f);CHKERRQ(ierr);
1747   PetscFunctionReturn(0);
1748 }
1749
1750 #undef __FUNCT__
1751 #define __FUNCT__ "DMDABCApplySymmetricCompression"
1752 static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff)
1753 {
1754   PetscErrorCode ierr;
1755   PetscInt       start,end,m;
1756   PetscInt       *unconstrained;
1757   PetscInt       cnt,i;
1758   Vec            x;
1759   PetscScalar    *_x;
1760   IS             is;
1761   VecScatter     scat;
1762
1763   PetscFunctionBeginUser;
1764   /* push bc's into f and A */
1765   ierr = VecDuplicate(f,&x);CHKERRQ(ierr);
1766   ierr = BCApply_EAST(elas_da,0,-1.0,A,x);CHKERRQ(ierr);
1767   ierr = BCApply_EAST(elas_da,1, 0.0,A,x);CHKERRQ(ierr);
1768   ierr = BCApply_WEST(elas_da,0,1.0,A,x);CHKERRQ(ierr);
1769   ierr = BCApply_WEST(elas_da,1,0.0,A,x);CHKERRQ(ierr);
1770
1771   /* define which dofs are not constrained */
1772   ierr = VecGetLocalSize(x,&m);CHKERRQ(ierr);
1773   ierr = PetscMalloc(sizeof(PetscInt)*m,&unconstrained);CHKERRQ(ierr);
1774   ierr = VecGetOwnershipRange(x,&start,&end);CHKERRQ(ierr);
1775   ierr = VecGetArray(x,&_x);CHKERRQ(ierr);
1776   cnt  = 0;
1777   for (i = 0; i < m; i++) {
1778     PetscReal val;
1779
1780     val = PetscRealPart(_x[i]);
1781     if (fabs(val) < 0.1) {
1782       unconstrained[cnt] = start + i;
1783       cnt++;
1784     }
1785   }
1786   ierr = VecRestoreArray(x,&_x);CHKERRQ(ierr);
1787
1788   ierr = ISCreateGeneral(PETSC_COMM_WORLD,cnt,unconstrained,PETSC_COPY_VALUES,&is);CHKERRQ(ierr);
1789   ierr = PetscFree(unconstrained);CHKERRQ(ierr);
1790
1791   /* define correction for dirichlet in the rhs */
1792   ierr = MatMult(A,x,f);CHKERRQ(ierr);
1793   ierr = VecScale(f,-1.0);CHKERRQ(ierr);
1794
1795   /* get new matrix */
1796   ierr = MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,AA);CHKERRQ(ierr);
1797   /* get new vector */
1798   ierr = MatGetVecs(*AA,NULL,ff);CHKERRQ(ierr);
1799
1800   ierr = VecScatterCreate(f,is,*ff,NULL,&scat);CHKERRQ(ierr);
1801   ierr = VecScatterBegin(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1802   ierr = VecScatterEnd(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1803
1804   {                             /* Constrain near-null space */
1805     PetscInt nvecs;
1806     const Vec *vecs;
1807     Vec *uvecs;
1808     PetscBool has_const;
1809     MatNullSpace mnull,unull;
1810     ierr = MatGetNearNullSpace(A,&mnull);CHKERRQ(ierr);
1811     ierr = MatNullSpaceGetVecs(mnull,&has_const,&nvecs,&vecs);CHKERRQ(ierr);
1812     ierr = VecDuplicateVecs(*ff,nvecs,&uvecs);CHKERRQ(ierr);
1813     for (i=0; i<nvecs; i++) {
1814       ierr = VecScatterBegin(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1815       ierr = VecScatterEnd(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1816     }
1817     ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)A),has_const,nvecs,uvecs,&unull);CHKERRQ(ierr);
1818     ierr = MatSetNearNullSpace(*AA,unull);CHKERRQ(ierr);
1819     ierr = MatNullSpaceDestroy(&unull);CHKERRQ(ierr);
1820     for (i=0; i<nvecs; i++) {
1821       ierr = VecDestroy(&uvecs[i]);CHKERRQ(ierr);
1822     }
1823     ierr = PetscFree(uvecs);CHKERRQ(ierr);
1824   }
1825
1826   ierr = VecScatterDestroy(&scat);CHKERRQ(ierr);
1827
1828   *dofs = is;
1829   ierr  = VecDestroy(&x);CHKERRQ(ierr);
1830   PetscFunctionReturn(0);
1831 }