code/ex49.c

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