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