Actual source code: ex43.c

petsc-3.5.1 2014-07-24
Report Typos and Errors
  1: static char help[] = "Solves the incompressible, variable viscosity stokes equation in 2d on the unit domain \n\
  2: using Q1Q1 elements, stabilized with Bochev's polynomial projection method. \n\
  3: The models defined utilise free slip boundary conditions on all sides. \n\
  4: Options: \n\
  5:      -mx : number elements in x-direciton \n\
  6:      -my : number elements in y-direciton \n\
  7:      -c_str : indicates the structure of the coefficients to use. \n\
  8:           -c_str 0 => Setup for an analytic solution with a vertical jump in viscosity. This problem is driven by the \n\
  9:                          forcing function f = (0, sin(n_z pi y)cos(pi x). \n\
 10:                          Parameters: \n\
 11:                               -solcx_eta0 : the viscosity to the left of the interface \n\
 12:                               -solcx_eta1 : the viscosity to the right of the interface \n\
 13:                               -solcx_xc : the location of the interface \n\
 14:                               -solcx_nz : the wavenumber in the y direction \n\
 15:           -c_str 1 => Setup for a rectangular blob located in the origin of the domain. \n\
 16:                          Parameters: \n\
 17:                               -sinker_eta0 : the viscosity of the background fluid \n\
 18:                               -sinker_eta1 : the viscosity of the blob \n\
 19:                               -sinker_dx : the width of the blob \n\
 20:                               -sinker_dy : the width of the blob \n\
 21:           -c_str 2 => Setup for a circular blob located in the origin of the domain. \n\
 22:                          Parameters: \n\
 23:                               -sinker_eta0 : the viscosity of the background fluid \n\
 24:                               -sinker_eta1 : the viscosity of the blob \n\
 25:                               -sinker_r : radius of the blob \n\
 26:           -c_str 3 => Circular and rectangular inclusion\n\
 27:                          Parameters as for cases 1 and 2 (above)\n\
 28:      -use_gp_coords : evaluate the viscosity and the body force at the global coordinates of the quadrature points.\n\
 29:      By default, eta and the body force are evaulated at the element center and applied as a constant over the entire element.\n";

 31: /* Contributed by Dave May */

 33: #include <petscksp.h>
 34: #include <petscdm.h>
 35: #include <petscdmda.h>

 37: /* A Maple-generated exact solution created by Mirko Velic (mirko.velic@sci.monash.edu.au) */
 38: #include "ex43-solcx.h"

 40: static PetscErrorCode DMDABCApplyFreeSlip(DM,Mat,Vec);


 43: #define NSD            2 /* number of spatial dimensions */
 44: #define NODES_PER_EL   4 /* nodes per element */
 45: #define U_DOFS         2 /* degrees of freedom per velocity node */
 46: #define P_DOFS         1 /* degrees of freedom per pressure node */
 47: #define GAUSS_POINTS   4

 49: /* cell based evaluation */
 50: typedef struct {
 51:   PetscScalar eta,fx,fy;
 52: } Coefficients;

 54: /* Gauss point based evaluation 8+4+4+4 = 20 */
 55: typedef struct {
 56:   PetscScalar gp_coords[2*GAUSS_POINTS];
 57:   PetscScalar eta[GAUSS_POINTS];
 58:   PetscScalar fx[GAUSS_POINTS];
 59:   PetscScalar fy[GAUSS_POINTS];
 60: } GaussPointCoefficients;

 62: typedef struct {
 63:   PetscScalar u_dof;
 64:   PetscScalar v_dof;
 65:   PetscScalar p_dof;
 66: } StokesDOF;


 69: /*

 71: D = [ 2.eta   0   0   ]
 72: [   0   2.eta 0   ]
 73: [   0     0   eta ]

 75: B = [ d_dx   0   ]
 76: [  0    d_dy ]
 77: [ d_dy  d_dx ]

 79: */

 81: /* FEM routines */
 82: /*
 83: Element: Local basis function ordering
 84: 1-----2
 85: |     |
 86: |     |
 87: 0-----3
 88: */
 89: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
 90: {
 91:   PetscScalar xi  = _xi[0];
 92:   PetscScalar eta = _xi[1];

 94:   Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
 95:   Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
 96:   Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
 97:   Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
 98: }

100: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
101: {
102:   PetscScalar xi  = _xi[0];
103:   PetscScalar eta = _xi[1];

105:   GNi[0][0] = -0.25*(1.0-eta);
106:   GNi[0][1] = -0.25*(1.0+eta);
107:   GNi[0][2] =   0.25*(1.0+eta);
108:   GNi[0][3] =   0.25*(1.0-eta);

110:   GNi[1][0] = -0.25*(1.0-xi);
111:   GNi[1][1] =   0.25*(1.0-xi);
112:   GNi[1][2] =   0.25*(1.0+xi);
113:   GNi[1][3] = -0.25*(1.0+xi);
114: }

116: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
117: {
118:   PetscScalar J00,J01,J10,J11,J;
119:   PetscScalar iJ00,iJ01,iJ10,iJ11;
120:   PetscInt    i;

122:   J00 = J01 = J10 = J11 = 0.0;
123:   for (i = 0; i < NODES_PER_EL; i++) {
124:     PetscScalar cx = coords[2*i+0];
125:     PetscScalar cy = coords[2*i+1];

127:     J00 = J00+GNi[0][i]*cx;      /* J_xx = dx/dxi */
128:     J01 = J01+GNi[0][i]*cy;      /* J_xy = dy/dxi */
129:     J10 = J10+GNi[1][i]*cx;      /* J_yx = dx/deta */
130:     J11 = J11+GNi[1][i]*cy;      /* J_yy = dy/deta */
131:   }
132:   J = (J00*J11)-(J01*J10);

134:   iJ00 =  J11/J;
135:   iJ01 = -J01/J;
136:   iJ10 = -J10/J;
137:   iJ11 =  J00/J;


140:   for (i = 0; i < NODES_PER_EL; i++) {
141:     GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
142:     GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
143:   }

145:   if (det_J != NULL) *det_J = J;
146: }

148: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
149: {
150:   *ngp         = 4;
151:   gp_xi[0][0]  = -0.57735026919;gp_xi[0][1] = -0.57735026919;
152:   gp_xi[1][0]  = -0.57735026919;gp_xi[1][1] =  0.57735026919;
153:   gp_xi[2][0]  =  0.57735026919;gp_xi[2][1] =  0.57735026919;
154:   gp_xi[3][0]  =  0.57735026919;gp_xi[3][1] = -0.57735026919;
155:   gp_weight[0] = 1.0;
156:   gp_weight[1] = 1.0;
157:   gp_weight[2] = 1.0;
158:   gp_weight[3] = 1.0;
159: }


162: /* procs to the left claim the ghost node as their element */
165: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
166: {
167:   PetscInt m,n,p,M,N,P;
168:   PetscInt sx,sy,sz;

171:   DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
172:   DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);

174:   if (mxl != NULL) {
175:     *mxl = m;
176:     if ((sx+m) == M) *mxl = m-1;  /* last proc */
177:   }
178:   if (myl != NULL) {
179:     *myl = n;
180:     if ((sy+n) == N) *myl = n-1;  /* last proc */
181:   }
182:   if (mzl != NULL) {
183:     *mzl = p;
184:     if ((sz+p) == P) *mzl = p-1;  /* last proc */
185:   }
186:   return(0);
187: }

191: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
192: {
193:   PetscInt si,sj,sk;

196:   DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);

198:   *sx = si;
199:   if (si) *sx = si+1;

201:   *sy = sj;
202:   if (sj) *sy = sj+1;

204:   if (sk) {
205:     *sz = sk;
206:     if (sk != 0) *sz = sk+1;
207:   }

209:   DMDAGetLocalElementSize(da,mx,my,mz);
210:   return(0);
211: }

213: /*
214: i,j are the element indices
215: The unknown is a vector quantity.
216: The s[].c is used to indicate the degree of freedom.
217: */
220: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
221: {
223:   /* velocity */
224:   /* node 0 */
225:   s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0;                         /* Vx0 */
226:   s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1;                         /* Vy0 */

228:   /* node 1 */
229:   s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0;                         /* Vx1 */
230:   s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1;                         /* Vy1 */

232:   /* node 2 */
233:   s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0;                         /* Vx2 */
234:   s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1;                         /* Vy2 */

236:   /* node 3 */
237:   s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0;                         /* Vx3 */
238:   s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1;                         /* Vy3 */


241:   /* pressure */
242:   s_p[0].i = i;s_p[0].j = j;s_p[0].c = 2;                         /* P0 */
243:   s_p[1].i = i;s_p[1].j = j+1;s_p[1].c = 2;                         /* P0 */
244:   s_p[2].i = i+1;s_p[2].j = j+1;s_p[2].c = 2;                         /* P1 */
245:   s_p[3].i = i+1;s_p[3].j = j;s_p[3].c = 2;                         /* P1 */
246:   return(0);
247: }

251: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
252: {
254:   PetscMPIInt    rank;
255:   PetscInt       proc_I,proc_J;
256:   PetscInt       cpu_x,cpu_y;
257:   PetscInt       local_mx,local_my;
258:   Vec            vlx,vly;
259:   PetscInt       *LX,*LY,i;
260:   PetscScalar    *_a;
261:   Vec            V_SEQ;
262:   VecScatter     ctx;

265:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

267:   DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);

269:   proc_J = rank/cpu_x;
270:   proc_I = rank-cpu_x*proc_J;

272:   PetscMalloc(sizeof(PetscInt)*cpu_x,&LX);
273:   PetscMalloc(sizeof(PetscInt)*cpu_y,&LY);

275:   DMDAGetLocalElementSize(da,&local_mx,&local_my,NULL);
276:   VecCreate(PETSC_COMM_WORLD,&vlx);
277:   VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
278:   VecSetFromOptions(vlx);

280:   VecCreate(PETSC_COMM_WORLD,&vly);
281:   VecSetSizes(vly,PETSC_DECIDE,cpu_y);
282:   VecSetFromOptions(vly);

284:   VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
285:   VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
286:   VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
287:   VecAssemblyBegin(vly);VecAssemblyEnd(vly);



291:   VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
292:   VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
293:   VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
294:   VecGetArray(V_SEQ,&_a);
295:   for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
296:   VecRestoreArray(V_SEQ,&_a);
297:   VecScatterDestroy(&ctx);
298:   VecDestroy(&V_SEQ);

300:   VecScatterCreateToAll(vly,&ctx,&V_SEQ);
301:   VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
302:   VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
303:   VecGetArray(V_SEQ,&_a);
304:   for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
305:   VecRestoreArray(V_SEQ,&_a);
306:   VecScatterDestroy(&ctx);
307:   VecDestroy(&V_SEQ);



311:   *_lx = LX;
312:   *_ly = LY;

314:   VecDestroy(&vlx);
315:   VecDestroy(&vly);
316:   return(0);
317: }

321: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
322: {
323:   DM             cda;
324:   Vec            coords;
325:   DMDACoor2d     **_coords;
326:   PetscInt       si,sj,nx,ny,i,j;
327:   FILE           *fp;
328:   char           fname[PETSC_MAX_PATH_LEN];
329:   PetscMPIInt    rank;

333:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
334:   PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
335:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
336:   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");

338:   PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);

340:   DMGetCoordinateDM(da,&cda);
341:   DMGetCoordinatesLocal(da,&coords);
342:   DMDAVecGetArray(cda,coords,&_coords);
343:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
344:   for (j = sj; j < sj+ny-1; j++) {
345:     for (i = si; i < si+nx-1; i++) {
346:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
347:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i].x),(double)PetscRealPart(_coords[j+1][i].y));
348:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i+1].x),(double)PetscRealPart(_coords[j+1][i+1].y));
349:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i+1].x),(double)PetscRealPart(_coords[j][i+1].y));
350:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
351:     }
352:   }
353:   DMDAVecRestoreArray(cda,coords,&_coords);

355:   PetscFClose(PETSC_COMM_SELF,fp);
356:   return(0);
357: }

361: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
362: {
363:   DM             cda;
364:   Vec            coords,local_fields;
365:   DMDACoor2d     **_coords;
366:   FILE           *fp;
367:   char           fname[PETSC_MAX_PATH_LEN];
368:   PetscMPIInt    rank;
369:   PetscInt       si,sj,nx,ny,i,j;
370:   PetscInt       n_dofs,d;
371:   PetscScalar    *_fields;

375:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
376:   PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
377:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
378:   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");

380:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
381:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
382:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
383:   for (d = 0; d < n_dofs; d++) {
384:     const char *field_name;
385:     DMDAGetFieldName(da,d,&field_name);
386:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
387:   }
388:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");


391:   DMGetCoordinateDM(da,&cda);
392:   DMGetCoordinatesLocal(da,&coords);
393:   DMDAVecGetArray(cda,coords,&_coords);
394:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

396:   DMCreateLocalVector(da,&local_fields);
397:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
398:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
399:   VecGetArray(local_fields,&_fields);


402:   for (j = sj; j < sj+ny; j++) {
403:     for (i = si; i < si+nx; i++) {
404:       PetscScalar coord_x,coord_y;
405:       PetscScalar field_d;

407:       coord_x = _coords[j][i].x;
408:       coord_y = _coords[j][i].y;

410:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
411:       for (d = 0; d < n_dofs; d++) {
412:         field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
413:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",(double)PetscRealPart(field_d));
414:       }
415:       PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
416:     }
417:   }
418:   VecRestoreArray(local_fields,&_fields);
419:   VecDestroy(&local_fields);

421:   DMDAVecRestoreArray(cda,coords,&_coords);

423:   PetscFClose(PETSC_COMM_SELF,fp);
424:   return(0);
425: }

429: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
430: {
431:   DM                     cda;
432:   Vec                    local_fields;
433:   FILE                   *fp;
434:   char                   fname[PETSC_MAX_PATH_LEN];
435:   PetscMPIInt            rank;
436:   PetscInt               si,sj,nx,ny,i,j,p;
437:   PetscInt               n_dofs,d;
438:   GaussPointCoefficients **_coefficients;
439:   PetscErrorCode         ierr;

442:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
443:   PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
444:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
445:   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");

447:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
448:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
449:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
450:   for (d = 0; d < n_dofs; d++) {
451:     const char *field_name;
452:     DMDAGetFieldName(da,d,&field_name);
453:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
454:   }
455:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");


458:   DMGetCoordinateDM(da,&cda);
459:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

461:   DMCreateLocalVector(da,&local_fields);
462:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
463:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
464:   DMDAVecGetArray(da,local_fields,&_coefficients);


467:   for (j = sj; j < sj+ny; j++) {
468:     for (i = si; i < si+nx; i++) {
469:       PetscScalar coord_x,coord_y;

471:       for (p = 0; p < GAUSS_POINTS; p++) {
472:         coord_x = _coefficients[j][i].gp_coords[2*p];
473:         coord_y = _coefficients[j][i].gp_coords[2*p+1];

475:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));

477:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e",(double)PetscRealPart(_coefficients[j][i].eta[p]),(double)PetscRealPart(_coefficients[j][i].fx[p]),(double)PetscRealPart(_coefficients[j][i].fy[p]));
478:         PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
479:       }
480:     }
481:   }
482:   DMDAVecRestoreArray(da,local_fields,&_coefficients);
483:   VecDestroy(&local_fields);

485:   PetscFClose(PETSC_COMM_SELF,fp);
486:   return(0);
487: }


490: static PetscInt ASS_MAP_wIwDI_uJuDJ(PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)
491: {
492:   PetscInt ij;
493:   PetscInt r,c,nc;

495:   nc = u_NPE*u_dof;

497:   r = w_dof*wi+wd;
498:   c = u_dof*ui+ud;

500:   ij = r*nc+c;

502:   return ij;
503: }

505: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
506: {
507:   PetscInt    ngp;
508:   PetscScalar gp_xi[GAUSS_POINTS][2];
509:   PetscScalar gp_weight[GAUSS_POINTS];
510:   PetscInt    p,i,j,k;
511:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
512:   PetscScalar J_p,tildeD[3];
513:   PetscScalar B[3][U_DOFS*NODES_PER_EL];


516:   /* define quadrature rule */
517:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

519:   /* evaluate integral */
520:   for (p = 0; p < ngp; p++) {
521:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
522:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);

524:     for (i = 0; i < NODES_PER_EL; i++) {
525:       PetscScalar d_dx_i = GNx_p[0][i];
526:       PetscScalar d_dy_i = GNx_p[1][i];

528:       B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
529:       B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
530:       B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
531:     }


534:     tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
535:     tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
536:     tildeD[2] =       gp_weight[p]*J_p*eta[p];

538:     /* form Bt tildeD B */
539:     /*
540:     Ke_ij = Bt_ik . D_kl . B_lj
541:     = B_ki . D_kl . B_lj
542:     = B_ki . D_kk . B_kj
543:     */
544:     for (i = 0; i < 8; i++) {
545:       for (j = 0; j < 8; j++) {
546:         for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
547:           Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
548:         }
549:       }
550:     }
551:   }
552: }

554: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
555: {
556:   PetscInt    ngp;
557:   PetscScalar gp_xi[GAUSS_POINTS][2];
558:   PetscScalar gp_weight[GAUSS_POINTS];
559:   PetscInt    p,i,j,di;
560:   PetscScalar Ni_p[NODES_PER_EL];
561:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
562:   PetscScalar J_p,fac;


565:   /* define quadrature rule */
566:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

568:   /* evaluate integral */
569:   for (p = 0; p < ngp; p++) {
570:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
571:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
572:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
573:     fac = gp_weight[p]*J_p;

575:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
576:       for (di = 0; di < NSD; di++) { /* u dofs */
577:         for (j = 0; j < 4; j++) {  /* p nodes, p dofs = 1 (ie no loop) */
578:           PetscInt IJ;
579:           /*     Ke[4*u_idx+j] = Ke[4*u_idx+j] - GNx_p[di][i] * Ni_p[j] * fac; */
580:           IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);

582:           Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
583:         }
584:       }
585:     }
586:   }
587: }

589: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
590: {
591:   PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
592:   PetscInt    i,j;
593:   PetscInt    nr_g,nc_g;

595:   PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
596:   FormGradientOperatorQ1(Ge,coords);

598:   nr_g = U_DOFS*NODES_PER_EL;
599:   nc_g = P_DOFS*NODES_PER_EL;

601:   for (i = 0; i < nr_g; i++) {
602:     for (j = 0; j < nc_g; j++) {
603:       De[nr_g*j+i] = Ge[nc_g*i+j];
604:     }
605:   }
606: }

608: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
609: {
610:   PetscInt    ngp;
611:   PetscScalar gp_xi[GAUSS_POINTS][2];
612:   PetscScalar gp_weight[GAUSS_POINTS];
613:   PetscInt    p,i,j;
614:   PetscScalar Ni_p[NODES_PER_EL];
615:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
616:   PetscScalar J_p,fac,eta_avg;


619:   /* define quadrature rule */
620:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

622:   /* evaluate integral */
623:   for (p = 0; p < ngp; p++) {
624:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
625:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
626:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
627:     fac = gp_weight[p]*J_p;

629:     for (i = 0; i < NODES_PER_EL; i++) {
630:       for (j = 0; j < NODES_PER_EL; j++) {
631:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
632:       }
633:     }
634:   }

636:   /* scale */
637:   eta_avg = 0.0;
638:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
639:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
640:   fac     = 1.0/eta_avg;
641:   for (i = 0; i < NODES_PER_EL; i++) {
642:     for (j = 0; j < NODES_PER_EL; j++) {
643:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
644:     }
645:   }
646: }

648: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
649: {
650:   PetscInt    ngp;
651:   PetscScalar gp_xi[GAUSS_POINTS][2];
652:   PetscScalar gp_weight[GAUSS_POINTS];
653:   PetscInt    p,i,j;
654:   PetscScalar Ni_p[NODES_PER_EL];
655:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
656:   PetscScalar J_p,fac,eta_avg;


659:   /* define quadrature rule */
660:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

662:   /* evaluate integral */
663:   for (p = 0; p < ngp; p++) {
664:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
665:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
666:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
667:     fac = gp_weight[p]*J_p;

669:     for (i = 0; i < NODES_PER_EL; i++) {
670:       for (j = 0; j < NODES_PER_EL; j++) {
671:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
672:       }
673:     }
674:   }

676:   /* scale */
677:   eta_avg = 0.0;
678:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
679:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
680:   fac     = 1.0/eta_avg;
681:   for (i = 0; i < NODES_PER_EL; i++) {
682:     for (j = 0; j < NODES_PER_EL; j++) {
683:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
684:     }
685:   }
686: }

688: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
689: {
690:   PetscInt    ngp;
691:   PetscScalar gp_xi[GAUSS_POINTS][2];
692:   PetscScalar gp_weight[GAUSS_POINTS];
693:   PetscInt    p,i;
694:   PetscScalar Ni_p[NODES_PER_EL];
695:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
696:   PetscScalar J_p,fac;


699:   /* define quadrature rule */
700:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

702:   /* evaluate integral */
703:   for (p = 0; p < ngp; p++) {
704:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
705:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
706:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
707:     fac = gp_weight[p]*J_p;

709:     for (i = 0; i < NODES_PER_EL; i++) {
710:       Fe[NSD*i]   += fac*Ni_p[i]*fx[p];
711:       Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
712:     }
713:   }
714: }

718: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
719: {
721:   /* get coords for the element */
722:   el_coords[NSD*0+0] = _coords[ej][ei].x;el_coords[NSD*0+1] = _coords[ej][ei].y;
723:   el_coords[NSD*1+0] = _coords[ej+1][ei].x;el_coords[NSD*1+1] = _coords[ej+1][ei].y;
724:   el_coords[NSD*2+0] = _coords[ej+1][ei+1].x;el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
725:   el_coords[NSD*3+0] = _coords[ej][ei+1].x;el_coords[NSD*3+1] = _coords[ej][ei+1].y;
726:   return(0);
727: }

731: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
732: {
733:   DM                     cda;
734:   Vec                    coords;
735:   DMDACoor2d             **_coords;
736:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
737:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
738:   PetscInt               sex,sey,mx,my;
739:   PetscInt               ei,ej;
740:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
741:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
742:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
743:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
744:   PetscScalar            el_coords[NODES_PER_EL*NSD];
745:   Vec                    local_properties;
746:   GaussPointCoefficients **props;
747:   PetscScalar            *prop_eta;
748:   PetscErrorCode         ierr;

751:   /* setup for coords */
752:   DMGetCoordinateDM(stokes_da,&cda);
753:   DMGetCoordinatesLocal(stokes_da,&coords);
754:   DMDAVecGetArray(cda,coords,&_coords);

756:   /* setup for coefficients */
757:   DMCreateLocalVector(properties_da,&local_properties);
758:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
759:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
760:   DMDAVecGetArray(properties_da,local_properties,&props);

762:   DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
763:   for (ej = sey; ej < sey+my; ej++) {
764:     for (ei = sex; ei < sex+mx; ei++) {
765:       /* get coords for the element */
766:       GetElementCoords(_coords,ei,ej,el_coords);

768:       /* get coefficients for the element */
769:       prop_eta = props[ej][ei].eta;

771:       /* initialise element stiffness matrix */
772:       PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
773:       PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
774:       PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
775:       PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);

777:       /* form element stiffness matrix */
778:       FormStressOperatorQ1(Ae,el_coords,prop_eta);
779:       FormGradientOperatorQ1(Ge,el_coords);
780:       FormDivergenceOperatorQ1(De,el_coords);
781:       FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);

783:       /* insert element matrix into global matrix */
784:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
785:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
786:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
787:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
788:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
789:     }
790:   }
791:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
792:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

794:   DMDAVecRestoreArray(cda,coords,&_coords);

796:   DMDAVecRestoreArray(properties_da,local_properties,&props);
797:   VecDestroy(&local_properties);
798:   return(0);
799: }

803: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
804: {
805:   DM                     cda;
806:   Vec                    coords;
807:   DMDACoor2d             **_coords;
808:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
809:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
810:   PetscInt               sex,sey,mx,my;
811:   PetscInt               ei,ej;
812:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
813:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
814:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
815:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
816:   PetscScalar            el_coords[NODES_PER_EL*NSD];
817:   Vec                    local_properties;
818:   GaussPointCoefficients **props;
819:   PetscScalar            *prop_eta;
820:   PetscErrorCode         ierr;

823:   /* setup for coords */
824:   DMGetCoordinateDM(stokes_da,&cda);
825:   DMGetCoordinatesLocal(stokes_da,&coords);
826:   DMDAVecGetArray(cda,coords,&_coords);

828:   /* setup for coefficients */
829:   DMCreateLocalVector(properties_da,&local_properties);
830:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
831:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
832:   DMDAVecGetArray(properties_da,local_properties,&props);

834:   DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
835:   for (ej = sey; ej < sey+my; ej++) {
836:     for (ei = sex; ei < sex+mx; ei++) {
837:       /* get coords for the element */
838:       GetElementCoords(_coords,ei,ej,el_coords);

840:       /* get coefficients for the element */
841:       prop_eta = props[ej][ei].eta;

843:       /* initialise element stiffness matrix */
844:       PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
845:       PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
846:       PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
847:       PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);


850:       /* form element stiffness matrix */
851:       FormStressOperatorQ1(Ae,el_coords,prop_eta);
852:       FormGradientOperatorQ1(Ge,el_coords);
853:       /*               FormDivergenceOperatorQ1(De, el_coords); */
854:       FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);

856:       /* insert element matrix into global matrix */
857:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
858:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
859:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
860:       /*     MatSetValuesStencil(A, NODES_PER_EL*P_DOFS,p_eqn, NODES_PER_EL*U_DOFS,u_eqn, De, ADD_VALUES); */
861:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
862:     }
863:   }
864:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
865:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

867:   DMDAVecRestoreArray(cda,coords,&_coords);

869:   DMDAVecRestoreArray(properties_da,local_properties,&props);
870:   VecDestroy(&local_properties);
871:   return(0);
872: }

876: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
877: {
878:   PetscInt n;

881:   for (n = 0; n < 4; n++) {
882:     fields_F[u_eqn[2*n].j][u_eqn[2*n].i].u_dof     = fields_F[u_eqn[2*n].j][u_eqn[2*n].i].u_dof+Fe_u[2*n];
883:     fields_F[u_eqn[2*n+1].j][u_eqn[2*n+1].i].v_dof = fields_F[u_eqn[2*n+1].j][u_eqn[2*n+1].i].v_dof+Fe_u[2*n+1];
884:     fields_F[p_eqn[n].j][p_eqn[n].i].p_dof         = fields_F[p_eqn[n].j][p_eqn[n].i].p_dof+Fe_p[n];
885:   }
886:   return(0);
887: }

891: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
892: {
893:   DM                     cda;
894:   Vec                    coords;
895:   DMDACoor2d             **_coords;
896:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
897:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
898:   PetscInt               sex,sey,mx,my;
899:   PetscInt               ei,ej;
900:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
901:   PetscScalar            He[NODES_PER_EL*P_DOFS];
902:   PetscScalar            el_coords[NODES_PER_EL*NSD];
903:   Vec                    local_properties;
904:   GaussPointCoefficients **props;
905:   PetscScalar            *prop_fx,*prop_fy;
906:   Vec                    local_F;
907:   StokesDOF              **ff;
908:   PetscErrorCode         ierr;

911:   /* setup for coords */
912:   DMGetCoordinateDM(stokes_da,&cda);
913:   DMGetCoordinatesLocal(stokes_da,&coords);
914:   DMDAVecGetArray(cda,coords,&_coords);

916:   /* setup for coefficients */
917:   DMGetLocalVector(properties_da,&local_properties);
918:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
919:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
920:   DMDAVecGetArray(properties_da,local_properties,&props);

922:   /* get acces to the vector */
923:   DMGetLocalVector(stokes_da,&local_F);
924:   VecZeroEntries(local_F);
925:   DMDAVecGetArray(stokes_da,local_F,&ff);


928:   DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
929:   for (ej = sey; ej < sey+my; ej++) {
930:     for (ei = sex; ei < sex+mx; ei++) {
931:       /* get coords for the element */
932:       GetElementCoords(_coords,ei,ej,el_coords);

934:       /* get coefficients for the element */
935:       prop_fx = props[ej][ei].fx;
936:       prop_fy = props[ej][ei].fy;

938:       /* initialise element stiffness matrix */
939:       PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
940:       PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);


943:       /* form element stiffness matrix */
944:       FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);

946:       /* insert element matrix into global matrix */
947:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);

949:       DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
950:     }
951:   }

953:   DMDAVecRestoreArray(stokes_da,local_F,&ff);
954:   DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
955:   DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
956:   DMRestoreLocalVector(stokes_da,&local_F);


959:   DMDAVecRestoreArray(cda,coords,&_coords);

961:   DMDAVecRestoreArray(properties_da,local_properties,&props);
962:   DMRestoreLocalVector(properties_da,&local_properties);
963:   return(0);
964: }

968: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,PetscInt mx,PetscInt my,DM *_da,Vec *_X)
969: {
970:   DM             da,cda;
971:   Vec            X,local_X;
972:   StokesDOF      **_stokes;
973:   Vec            coords;
974:   DMDACoor2d     **_coords;
975:   PetscInt       si,sj,ei,ej,i,j;

979:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
980:                       mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,&da);
981:   DMDASetFieldName(da,0,"anlytic_Vx");
982:   DMDASetFieldName(da,1,"anlytic_Vy");
983:   DMDASetFieldName(da,2,"analytic_P");


986:   DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.,0.);


989:   DMGetCoordinatesLocal(da,&coords);
990:   DMGetCoordinateDM(da,&cda);
991:   DMDAVecGetArray(cda,coords,&_coords);

993:   DMCreateGlobalVector(da,&X);
994:   DMCreateLocalVector(da,&local_X);
995:   DMDAVecGetArray(da,local_X,&_stokes);

997:   DMDAGetGhostCorners(da,&si,&sj,0,&ei,&ej,0);
998:   for (j = sj; j < sj+ej; j++) {
999:     for (i = si; i < si+ei; i++) {
1000:       PetscReal pos[2],pressure,vel[2],total_stress[3],strain_rate[3];

1002:       pos[0] = PetscRealPart(_coords[j][i].x);
1003:       pos[1] = PetscRealPart(_coords[j][i].y);

1005:       evaluate_solCx(pos,eta0,eta1,xc,nz,vel,&pressure,total_stress,strain_rate);

1007:       _stokes[j][i].u_dof = vel[0];
1008:       _stokes[j][i].v_dof = vel[1];
1009:       _stokes[j][i].p_dof = pressure;
1010:     }
1011:   }
1012:   DMDAVecRestoreArray(da,local_X,&_stokes);
1013:   DMDAVecRestoreArray(cda,coords,&_coords);

1015:   DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1016:   DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);

1018:   VecDestroy(&local_X);

1020:   *_da = da;
1021:   *_X  = X;
1022:   return(0);
1023: }

1027: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
1028: {
1030:   /* get the nodal fields */
1031:   nodal_fields[0].u_dof = fields[ej][ei].u_dof;nodal_fields[0].v_dof = fields[ej][ei].v_dof;nodal_fields[0].p_dof = fields[ej][ei].p_dof;
1032:   nodal_fields[1].u_dof = fields[ej+1][ei].u_dof;nodal_fields[1].v_dof = fields[ej+1][ei].v_dof;nodal_fields[1].p_dof = fields[ej+1][ei].p_dof;
1033:   nodal_fields[2].u_dof = fields[ej+1][ei+1].u_dof;nodal_fields[2].v_dof = fields[ej+1][ei+1].v_dof;nodal_fields[2].p_dof = fields[ej+1][ei+1].p_dof;
1034:   nodal_fields[3].u_dof = fields[ej][ei+1].u_dof;nodal_fields[3].v_dof = fields[ej][ei+1].v_dof;nodal_fields[3].p_dof = fields[ej][ei+1].p_dof;
1035:   return(0);
1036: }

1040: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
1041: {
1042:   DM          cda;
1043:   Vec         coords,X_analytic_local,X_local;
1044:   DMDACoor2d  **_coords;
1045:   PetscInt    sex,sey,mx,my;
1046:   PetscInt    ei,ej;
1047:   PetscScalar el_coords[NODES_PER_EL*NSD];
1048:   StokesDOF   **stokes_analytic,**stokes;
1049:   StokesDOF   stokes_analytic_e[4],stokes_e[4];

1051:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1052:   PetscScalar    Ni_p[NODES_PER_EL];
1053:   PetscInt       ngp;
1054:   PetscScalar    gp_xi[GAUSS_POINTS][2];
1055:   PetscScalar    gp_weight[GAUSS_POINTS];
1056:   PetscInt       p,i;
1057:   PetscScalar    J_p,fac;
1058:   PetscScalar    h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1059:   PetscInt       M;
1060:   PetscReal      xymin[2],xymax[2];

1064:   /* define quadrature rule */
1065:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

1067:   /* setup for coords */
1068:   DMGetCoordinateDM(stokes_da,&cda);
1069:   DMGetCoordinatesLocal(stokes_da,&coords);
1070:   DMDAVecGetArray(cda,coords,&_coords);

1072:   /* setup for analytic */
1073:   DMCreateLocalVector(stokes_da,&X_analytic_local);
1074:   DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1075:   DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1076:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1078:   /* setup for solution */
1079:   DMCreateLocalVector(stokes_da,&X_local);
1080:   DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1081:   DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1082:   DMDAVecGetArray(stokes_da,X_local,&stokes);

1084:   DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1085:   DMDAGetBoundingBox(stokes_da,xymin,xymax);

1087:   h = (xymax[0]-xymin[0])/((PetscReal)M);

1089:   tp_L2 = tu_L2 = tu_H1 = 0.0;

1091:   DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
1092:   for (ej = sey; ej < sey+my; ej++) {
1093:     for (ei = sex; ei < sex+mx; ei++) {
1094:       /* get coords for the element */
1095:       GetElementCoords(_coords,ei,ej,el_coords);
1096:       StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1097:       StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);

1099:       /* evaluate integral */
1100:       p_e_L2 = 0.0;
1101:       u_e_L2 = 0.0;
1102:       u_e_H1 = 0.0;
1103:       for (p = 0; p < ngp; p++) {
1104:         ConstructQ12D_Ni(gp_xi[p],Ni_p);
1105:         ConstructQ12D_GNi(gp_xi[p],GNi_p);
1106:         ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1107:         fac = gp_weight[p]*J_p;

1109:         for (i = 0; i < NODES_PER_EL; i++) {
1110:           PetscScalar u_error,v_error;

1112:           p_e_L2 = p_e_L2+fac*Ni_p[i]*(stokes_e[i].p_dof-stokes_analytic_e[i].p_dof)*(stokes_e[i].p_dof-stokes_analytic_e[i].p_dof);

1114:           u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1115:           v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1116:           u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);

1118:           u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error              /* du/dx */
1119:                                +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error               /* du/dy */
1120:                                +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error               /* dv/dx */
1121:                                +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error);             /* dv/dy */
1122:         }
1123:       }

1125:       tp_L2 += p_e_L2;
1126:       tu_L2 += u_e_L2;
1127:       tu_H1 += u_e_H1;
1128:     }
1129:   }
1130:   MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1131:   MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1132:   MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1133:   p_L2 = PetscSqrtScalar(p_L2);
1134:   u_L2 = PetscSqrtScalar(u_L2);
1135:   u_H1 = PetscSqrtScalar(u_H1);

1137:   PetscPrintf(PETSC_COMM_WORLD,"%1.4e   %1.4e   %1.4e   %1.4e \n",(double)PetscRealPart(h),(double)PetscRealPart(p_L2),(double)PetscRealPart(u_L2),(double)PetscRealPart(u_H1));


1140:   DMDAVecRestoreArray(cda,coords,&_coords);

1142:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1143:   VecDestroy(&X_analytic_local);
1144:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1145:   VecDestroy(&X_local);
1146:   return(0);
1147: }

1151: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1152: {
1153:   DM                     da_Stokes,da_prop;
1154:   PetscInt               u_dof,p_dof,dof,stencil_width;
1155:   Mat                    A,B;
1156:   PetscInt               mxl,myl;
1157:   DM                     prop_cda,vel_cda;
1158:   Vec                    prop_coords,vel_coords;
1159:   PetscInt               si,sj,nx,ny,i,j,p;
1160:   Vec                    f,X;
1161:   PetscInt               prop_dof,prop_stencil_width;
1162:   Vec                    properties,l_properties;
1163:   PetscReal              dx,dy;
1164:   PetscInt               M,N;
1165:   DMDACoor2d             **_prop_coords,**_vel_coords;
1166:   GaussPointCoefficients **element_props;
1167:   PetscInt               its;
1168:   KSP                    ksp_S;
1169:   PetscInt               coefficient_structure = 0;
1170:   PetscInt               cpu_x,cpu_y,*lx = NULL,*ly = NULL;
1171:   PetscBool              use_gp_coords = PETSC_FALSE,set;
1172:   char                   filename[PETSC_MAX_PATH_LEN];
1173:   PetscErrorCode         ierr;

1176:   /* Generate the da for velocity and pressure */
1177:   /*
1178:   We use Q1 elements for the temperature.
1179:   FEM has a 9-point stencil (BOX) or connectivity pattern
1180:   Num nodes in each direction is mx+1, my+1
1181:   */
1182:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1183:   p_dof         = P_DOFS; /* p - pressure */
1184:   dof           = u_dof+p_dof;
1185:   stencil_width = 1;
1186:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1187:                                mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&da_Stokes);
1188:   DMDASetFieldName(da_Stokes,0,"Vx");
1189:   DMDASetFieldName(da_Stokes,1,"Vy");
1190:   DMDASetFieldName(da_Stokes,2,"P");

1192:   /* unit box [0,1] x [0,1] */
1193:   DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.,0.);


1196:   /* Generate element properties, we will assume all material properties are constant over the element */
1197:   /* local number of elements */
1198:   DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,NULL);

1200:   /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!!  */
1201:   DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
1202:   DMDAGetElementOwnershipRanges2d(da_Stokes,&lx,&ly);

1204:   prop_dof           = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1205:   prop_stencil_width = 0;
1206:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1207:                                     mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
1208:   PetscFree(lx);
1209:   PetscFree(ly);

1211:   /* define centroid positions */
1212:   DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1213:   dx   = 1.0/((PetscReal)(M));
1214:   dy   = 1.0/((PetscReal)(N));

1216:   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);

1218:   /* define coefficients */
1219:   PetscOptionsGetInt(NULL,"-c_str",&coefficient_structure,NULL);
1220:   /*     PetscPrintf(PETSC_COMM_WORLD, "Using coeficient structure %D \n", coefficient_structure); */

1222:   DMCreateGlobalVector(da_prop,&properties);
1223:   DMCreateLocalVector(da_prop,&l_properties);
1224:   DMDAVecGetArray(da_prop,l_properties,&element_props);

1226:   DMGetCoordinateDM(da_prop,&prop_cda);
1227:   DMGetCoordinatesLocal(da_prop,&prop_coords);
1228:   DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);

1230:   DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);

1232:   DMGetCoordinateDM(da_Stokes,&vel_cda);
1233:   DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1234:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);


1237:   /* interpolate the coordinates */
1238:   for (j = sj; j < sj+ny; j++) {
1239:     for (i = si; i < si+nx; i++) {
1240:       PetscInt    ngp;
1241:       PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1242:       PetscScalar el_coords[8];

1244:       GetElementCoords(_vel_coords,i,j,el_coords);
1245:       ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

1247:       for (p = 0; p < GAUSS_POINTS; p++) {
1248:         PetscScalar gp_x,gp_y;
1249:         PetscInt    n;
1250:         PetscScalar xi_p[2],Ni_p[4];

1252:         xi_p[0] = gp_xi[p][0];
1253:         xi_p[1] = gp_xi[p][1];
1254:         ConstructQ12D_Ni(xi_p,Ni_p);

1256:         gp_x = 0.0;
1257:         gp_y = 0.0;
1258:         for (n = 0; n < NODES_PER_EL; n++) {
1259:           gp_x = gp_x+Ni_p[n]*el_coords[2*n];
1260:           gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1261:         }
1262:         element_props[j][i].gp_coords[2*p]   = gp_x;
1263:         element_props[j][i].gp_coords[2*p+1] = gp_y;
1264:       }
1265:     }
1266:   }

1268:   /* define the coefficients */
1269:   PetscOptionsGetBool(NULL,"-use_gp_coords",&use_gp_coords,0);

1271:   for (j = sj; j < sj+ny; j++) {
1272:     for (i = si; i < si+nx; i++) {
1273:       PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1274:       PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1275:       PetscReal coord_x,coord_y;

1277:       if (coefficient_structure == 0) {
1278:         PetscReal opts_eta0,opts_eta1,opts_xc;
1279:         PetscInt  opts_nz;

1281:         opts_eta0 = 1.0;
1282:         opts_eta1 = 1.0;
1283:         opts_xc   = 0.5;
1284:         opts_nz   = 1;

1286:         PetscOptionsGetReal(NULL,"-solcx_eta0",&opts_eta0,0);
1287:         PetscOptionsGetReal(NULL,"-solcx_eta1",&opts_eta1,0);
1288:         PetscOptionsGetReal(NULL,"-solcx_xc",&opts_xc,0);
1289:         PetscOptionsGetInt(NULL,"-solcx_nz",&opts_nz,0);

1291:         for (p = 0; p < GAUSS_POINTS; p++) {
1292:           coord_x = centroid_x;
1293:           coord_y = centroid_y;
1294:           if (use_gp_coords) {
1295:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1296:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1297:           }


1300:           element_props[j][i].eta[p] = opts_eta0;
1301:           if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;

1303:           element_props[j][i].fx[p] = 0.0;
1304:           element_props[j][i].fy[p] = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1305:         }
1306:       } else if (coefficient_structure == 1) { /* square sinker */
1307:         PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;

1309:         opts_eta0 = 1.0;
1310:         opts_eta1 = 1.0;
1311:         opts_dx   = 0.50;
1312:         opts_dy   = 0.50;

1314:         PetscOptionsGetReal(NULL,"-sinker_eta0",&opts_eta0,0);
1315:         PetscOptionsGetReal(NULL,"-sinker_eta1",&opts_eta1,0);
1316:         PetscOptionsGetReal(NULL,"-sinker_dx",&opts_dx,0);
1317:         PetscOptionsGetReal(NULL,"-sinker_dy",&opts_dy,0);


1320:         for (p = 0; p < GAUSS_POINTS; p++) {
1321:           coord_x = centroid_x;
1322:           coord_y = centroid_y;
1323:           if (use_gp_coords) {
1324:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1325:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1326:           }

1328:           element_props[j][i].eta[p] = opts_eta0;
1329:           element_props[j][i].fx[p]  = 0.0;
1330:           element_props[j][i].fy[p]  = 0.0;

1332:           if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1333:             if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1334:               element_props[j][i].eta[p] =  opts_eta1;
1335:               element_props[j][i].fx[p]  =  0.0;
1336:               element_props[j][i].fy[p]  = -1.0;
1337:             }
1338:           }
1339:         }
1340:       } else if (coefficient_structure == 2) { /* circular sinker */
1341:         PetscReal opts_eta0,opts_eta1,opts_r,radius2;

1343:         opts_eta0 = 1.0;
1344:         opts_eta1 = 1.0;
1345:         opts_r    = 0.25;

1347:         PetscOptionsGetReal(NULL,"-sinker_eta0",&opts_eta0,0);
1348:         PetscOptionsGetReal(NULL,"-sinker_eta1",&opts_eta1,0);
1349:         PetscOptionsGetReal(NULL,"-sinker_r",&opts_r,0);

1351:         for (p = 0; p < GAUSS_POINTS; p++) {
1352:           coord_x = centroid_x;
1353:           coord_y = centroid_y;
1354:           if (use_gp_coords) {
1355:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1356:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1357:           }

1359:           element_props[j][i].eta[p] = opts_eta0;
1360:           element_props[j][i].fx[p]  = 0.0;
1361:           element_props[j][i].fy[p]  = 0.0;

1363:           radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1364:           if (radius2 < opts_r*opts_r) {
1365:             element_props[j][i].eta[p] =  opts_eta1;
1366:             element_props[j][i].fx[p]  =  0.0;
1367:             element_props[j][i].fy[p]  = -1.0;
1368:           }
1369:         }
1370:       } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1371:         PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;

1373:         opts_eta0 = 1.0;
1374:         opts_eta1 = 1.0;
1375:         opts_r    = 0.25;
1376:         opts_c0x  = 0.35;       /* circle center */
1377:         opts_c0y  = 0.35;
1378:         opts_s0x  = 0.7;       /* square center */
1379:         opts_s0y  = 0.7;
1380:         opts_dx   = 0.25;
1381:         opts_dy   = 0.25;
1382:         opts_phi  = 25;

1384:         PetscOptionsGetReal(NULL,"-sinker_eta0",&opts_eta0,0);
1385:         PetscOptionsGetReal(NULL,"-sinker_eta1",&opts_eta1,0);
1386:         PetscOptionsGetReal(NULL,"-sinker_r",&opts_r,0);
1387:         PetscOptionsGetReal(NULL,"-sinker_c0x",&opts_c0x,0);
1388:         PetscOptionsGetReal(NULL,"-sinker_c0y",&opts_c0y,0);
1389:         PetscOptionsGetReal(NULL,"-sinker_s0x",&opts_s0x,0);
1390:         PetscOptionsGetReal(NULL,"-sinker_s0y",&opts_s0y,0);
1391:         PetscOptionsGetReal(NULL,"-sinker_dx",&opts_dx,0);
1392:         PetscOptionsGetReal(NULL,"-sinker_dy",&opts_dy,0);
1393:         PetscOptionsGetReal(NULL,"-sinker_phi",&opts_phi,0);
1394:         opts_phi *= PETSC_PI / 180;

1396:         for (p = 0; p < GAUSS_POINTS; p++) {
1397:           coord_x = centroid_x;
1398:           coord_y = centroid_y;
1399:           if (use_gp_coords) {
1400:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1401:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1402:           }

1404:           element_props[j][i].eta[p] = opts_eta0;
1405:           element_props[j][i].fx[p]  = 0.0;
1406:           element_props[j][i].fy[p]  = -0.2;

1408:           radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1409:           if (radius2 < opts_r*opts_r
1410:               || (PetscAbs(+(coord_x - opts_s0x)*PetscCosReal(opts_phi) + (coord_y - opts_s0y)*PetscSinReal(opts_phi)) < opts_dx/2
1411:                   && PetscAbs(-(coord_x - opts_s0x)*PetscSinReal(opts_phi) + (coord_y - opts_s0y)*PetscCosReal(opts_phi)) < opts_dy/2)) {
1412:             element_props[j][i].eta[p] =  opts_eta1;
1413:             element_props[j][i].fx[p]  =  0.0;
1414:             element_props[j][i].fy[p]  = -1.0;
1415:           }
1416:         }
1417:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1418:     }
1419:   }
1420:   DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);

1422:   DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);

1424:   DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1425:   DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1426:   DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);


1429:   DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1430:   DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for Stokes eqn.","properties");


1433:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1434:   DMSetMatType(da_Stokes,MATAIJ);
1435:   DMCreateMatrix(da_Stokes,&A);
1436:   DMCreateMatrix(da_Stokes,&B);
1437:   DMCreateGlobalVector(da_Stokes,&f);
1438:   DMCreateGlobalVector(da_Stokes,&X);

1440:   /* assemble A11 */
1441:   MatZeroEntries(A);
1442:   MatZeroEntries(B);
1443:   VecZeroEntries(f);

1445:   AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1446:   AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1447:   /* build force vector */
1448:   AssembleF_Stokes(f,da_Stokes,da_prop,properties);

1450:   DMDABCApplyFreeSlip(da_Stokes,A,f);
1451:   DMDABCApplyFreeSlip(da_Stokes,B,NULL);

1453:   /* SOLVE */
1454:   KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1455:   KSPSetOptionsPrefix(ksp_S,"stokes_");
1456:   KSPSetOperators(ksp_S,A,B);
1457:   KSPSetDM(ksp_S,da_Stokes);
1458:   KSPSetDMActive(ksp_S,PETSC_FALSE);
1459:   KSPSetFromOptions(ksp_S);
1460:   {
1461:     PC             pc;
1462:     const PetscInt ufields[] = {0,1},pfields[1] = {2};
1463:     KSPGetPC(ksp_S,&pc);
1464:     PCFieldSplitSetBlockSize(pc,3);
1465:     PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1466:     PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1467:   }

1469:   KSPSolve(ksp_S,f,X);

1471:   PetscOptionsGetString(NULL,"-o",filename,sizeof(filename),&set);
1472:   if (set) {
1473:     char        *ext;
1474:     PetscViewer viewer;
1475:     PetscBool   flg;
1476:     PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1477:     PetscStrrchr(filename,'.',&ext);
1478:     PetscStrcmp("vts",ext,&flg);
1479:     if (flg) {
1480:       PetscViewerSetType(viewer,PETSCVIEWERVTK);
1481:     } else {
1482:       PetscViewerSetType(viewer,PETSCVIEWERBINARY);
1483:     }
1484:     PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1485:     PetscViewerFileSetName(viewer,filename);
1486:     VecView(X,viewer);
1487:     PetscViewerDestroy(&viewer);
1488:   }
1489:   DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");

1491:   KSPGetIterationNumber(ksp_S,&its);

1493:   if (coefficient_structure == 0) {
1494:     PetscReal opts_eta0,opts_eta1,opts_xc;
1495:     PetscInt  opts_nz,N;
1496:     DM        da_Stokes_analytic;
1497:     Vec       X_analytic;
1498:     PetscReal nrm1[3],nrm2[3],nrmI[3];

1500:     opts_eta0 = 1.0;
1501:     opts_eta1 = 1.0;
1502:     opts_xc   = 0.5;
1503:     opts_nz   = 1;

1505:     PetscOptionsGetReal(NULL,"-solcx_eta0",&opts_eta0,0);
1506:     PetscOptionsGetReal(NULL,"-solcx_eta1",&opts_eta1,0);
1507:     PetscOptionsGetReal(NULL,"-solcx_xc",&opts_xc,0);
1508:     PetscOptionsGetInt(NULL,"-solcx_nz",&opts_nz,0);


1511:     DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1512:     DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");

1514:     DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);


1517:     VecAXPY(X_analytic,-1.0,X);
1518:     VecGetSize(X_analytic,&N);
1519:     N    = N/3;

1521:     VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1522:     VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1523:     VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);

1525:     VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1526:     VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1527:     VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);

1529:     VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1530:     VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1531:     VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);

1533:     DMDestroy(&da_Stokes_analytic);
1534:     VecDestroy(&X_analytic);
1535:   }


1538:   KSPDestroy(&ksp_S);
1539:   VecDestroy(&X);
1540:   VecDestroy(&f);
1541:   MatDestroy(&A);
1542:   MatDestroy(&B);

1544:   DMDestroy(&da_Stokes);
1545:   DMDestroy(&da_prop);

1547:   VecDestroy(&properties);
1548:   VecDestroy(&l_properties);
1549:   return(0);
1550: }

1554: int main(int argc,char **args)
1555: {
1557:   PetscInt       mx,my;

1559:   PetscInitialize(&argc,&args,(char*)0,help);

1561:   mx   = my = 10;
1562:   PetscOptionsGetInt(NULL,"-mx",&mx,NULL);
1563:   PetscOptionsGetInt(NULL,"-my",&my,NULL);

1565:   solve_stokes_2d_coupled(mx,my);

1567:   PetscFinalize();
1568:   return 0;
1569: }

1571: /* -------------------------- helpers for boundary conditions -------------------------------- */

1575: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1576: {
1577:   DM                     cda;
1578:   Vec                    coords;
1579:   PetscInt               si,sj,nx,ny,i,j;
1580:   PetscInt               M,N;
1581:   DMDACoor2d             **_coords;
1582:   const PetscInt         *g_idx;
1583:   PetscInt               *bc_global_ids;
1584:   PetscScalar            *bc_vals;
1585:   PetscInt               nbcs;
1586:   PetscInt               n_dofs;
1587:   PetscErrorCode         ierr;
1588:   ISLocalToGlobalMapping ltogm;

1591:   /* enforce bc's */
1592:   DMGetLocalToGlobalMapping(da,&ltogm);
1593:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1595:   DMGetCoordinateDM(da,&cda);
1596:   DMGetCoordinatesLocal(da,&coords);
1597:   DMDAVecGetArray(cda,coords,&_coords);
1598:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1599:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1601:   /* --- */

1603:   PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1604:   PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);

1606:   /* init the entries to -1 so VecSetValues will ignore them */
1607:   for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;

1609:   i = nx-1;
1610:   for (j = 0; j < ny; j++) {
1611:     PetscInt local_id;

1613:     local_id = i+j*nx;

1615:     bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];

1617:     bc_vals[j] =  bc_val;
1618:   }
1619:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1620:   nbcs = 0;
1621:   if ((si+nx) == (M)) nbcs = ny;

1623:   if (b != NULL) {
1624:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1625:     VecAssemblyBegin(b);
1626:     VecAssemblyEnd(b);
1627:   }
1628:   if (A != NULL) {
1629:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1630:   }


1633:   PetscFree(bc_vals);
1634:   PetscFree(bc_global_ids);

1636:   DMDAVecRestoreArray(cda,coords,&_coords);
1637:   return(0);
1638: }

1642: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1643: {
1644:   DM                     cda;
1645:   Vec                    coords;
1646:   PetscInt               si,sj,nx,ny,i,j;
1647:   PetscInt               M,N;
1648:   DMDACoor2d             **_coords;
1649:   const PetscInt         *g_idx;
1650:   PetscInt               *bc_global_ids;
1651:   PetscScalar            *bc_vals;
1652:   PetscInt               nbcs;
1653:   PetscInt               n_dofs;
1654:   PetscErrorCode         ierr;
1655:   ISLocalToGlobalMapping ltogm;

1658:   /* enforce bc's */
1659:   DMGetLocalToGlobalMapping(da,&ltogm);
1660:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1662:   DMGetCoordinateDM(da,&cda);
1663:   DMGetCoordinatesLocal(da,&coords);
1664:   DMDAVecGetArray(cda,coords,&_coords);
1665:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1666:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1668:   /* --- */

1670:   PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1671:   PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);

1673:   /* init the entries to -1 so VecSetValues will ignore them */
1674:   for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;

1676:   i = 0;
1677:   for (j = 0; j < ny; j++) {
1678:     PetscInt local_id;

1680:     local_id = i+j*nx;

1682:     bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];

1684:     bc_vals[j] =  bc_val;
1685:   }
1686:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1687:   nbcs = 0;
1688:   if (si == 0) nbcs = ny;

1690:   if (b != NULL) {
1691:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1692:     VecAssemblyBegin(b);
1693:     VecAssemblyEnd(b);
1694:   }
1695:   if (A != NULL) {
1696:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1697:   }


1700:   PetscFree(bc_vals);
1701:   PetscFree(bc_global_ids);

1703:   DMDAVecRestoreArray(cda,coords,&_coords);
1704:   return(0);
1705: }

1709: static PetscErrorCode BCApply_NORTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1710: {
1711:   DM                     cda;
1712:   Vec                    coords;
1713:   PetscInt               si,sj,nx,ny,i,j;
1714:   PetscInt               M,N;
1715:   DMDACoor2d             **_coords;
1716:   const PetscInt         *g_idx;
1717:   PetscInt               *bc_global_ids;
1718:   PetscScalar            *bc_vals;
1719:   PetscInt               nbcs;
1720:   PetscInt               n_dofs;
1721:   PetscErrorCode         ierr;
1722:   ISLocalToGlobalMapping ltogm;

1725:   /* enforce bc's */
1726:   DMGetLocalToGlobalMapping(da,&ltogm);
1727:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1729:   DMGetCoordinateDM(da,&cda);
1730:   DMGetCoordinatesLocal(da,&coords);
1731:   DMDAVecGetArray(cda,coords,&_coords);
1732:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1733:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1735:   /* --- */

1737:   PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1738:   PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);

1740:   /* init the entries to -1 so VecSetValues will ignore them */
1741:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

1743:   j = ny-1;
1744:   for (i = 0; i < nx; i++) {
1745:     PetscInt local_id;

1747:     local_id = i+j*nx;

1749:     bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];

1751:     bc_vals[i] =  bc_val;
1752:   }
1753:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1754:   nbcs = 0;
1755:   if ((sj+ny) == (N)) nbcs = nx;

1757:   if (b != NULL) {
1758:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1759:     VecAssemblyBegin(b);
1760:     VecAssemblyEnd(b);
1761:   }
1762:   if (A != NULL) {
1763:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1764:   }


1767:   PetscFree(bc_vals);
1768:   PetscFree(bc_global_ids);

1770:   DMDAVecRestoreArray(cda,coords,&_coords);
1771:   return(0);
1772: }

1776: static PetscErrorCode BCApply_SOUTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1777: {
1778:   DM                     cda;
1779:   Vec                    coords;
1780:   PetscInt               si,sj,nx,ny,i,j;
1781:   PetscInt               M,N;
1782:   DMDACoor2d             **_coords;
1783:   const PetscInt         *g_idx;
1784:   PetscInt               *bc_global_ids;
1785:   PetscScalar            *bc_vals;
1786:   PetscInt               nbcs;
1787:   PetscInt               n_dofs;
1788:   PetscErrorCode         ierr;
1789:   ISLocalToGlobalMapping ltogm;

1792:   /* enforce bc's */
1793:   DMGetLocalToGlobalMapping(da,&ltogm);
1794:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1796:   DMGetCoordinateDM(da,&cda);
1797:   DMGetCoordinatesLocal(da,&coords);
1798:   DMDAVecGetArray(cda,coords,&_coords);
1799:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1800:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1802:   /* --- */

1804:   PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1805:   PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);

1807:   /* init the entries to -1 so VecSetValues will ignore them */
1808:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

1810:   j = 0;
1811:   for (i = 0; i < nx; i++) {
1812:     PetscInt local_id;

1814:     local_id = i+j*nx;

1816:     bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];

1818:     bc_vals[i] =  bc_val;
1819:   }
1820:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1821:   nbcs = 0;
1822:   if (sj == 0) nbcs = nx;

1824:   if (b != NULL) {
1825:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1826:     VecAssemblyBegin(b);
1827:     VecAssemblyEnd(b);
1828:   }
1829:   if (A != NULL) {
1830:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1831:   }


1834:   PetscFree(bc_vals);
1835:   PetscFree(bc_global_ids);

1837:   DMDAVecRestoreArray(cda,coords,&_coords);
1838:   return(0);
1839: }

1843: /*
1844: Free slip sides.
1845: */
1848: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1849: {

1853:   BCApply_NORTH(da_Stokes,1,0.0,A,f);
1854:   BCApply_EAST(da_Stokes,0,0.0,A,f);
1855:   BCApply_SOUTH(da_Stokes,1,0.0,A,f);
1856:   BCApply_WEST(da_Stokes,0,0.0,A,f);
1857:   return(0);
1858: }