Actual source code: ex42.c

petsc-3.3-p7 2013-05-11
  1: static char help[] =
  2:   "Solves the incompressible, variable viscosity stokes equation in 3d using Q1Q1 elements, \n\
  3: stabilized with Bochev's polynomial projection method. Note that implementation here assumes \n\
  4: all boundaries are free-slip, i.e. zero normal flow and zero tangential stress \n\
  5:      -mx : number elements in x-direction \n\
  6:      -my : number elements in y-direction \n\
  7:      -mz : number elements in z-direction \n\
  8:      -stokes_ksp_monitor_blocks : active monitor for each component u,v,w,p \n\
  9:      -model : defines viscosity and forcing function. Choose either: 0 (isoviscous), 1 (manufactured-broken), 2 (solcx-2d), 3 (sinker) \n";

 11: /* Contributed by Dave May */

 13:  #include petscksp.h
 14:  #include petscpcmg.h
 15:  #include petscdmda.h

 17: #define PROFILE_TIMING
 18: #define ASSEMBLE_LOWER_TRIANGULAR

 20: #define NSD            3 /* number of spatial dimensions */
 21: #define NODES_PER_EL   8 /* nodes per element */
 22: #define U_DOFS         3 /* degrees of freedom per velocity node */
 23: #define P_DOFS         1 /* degrees of freedom per pressure node */
 24: #define GAUSS_POINTS   8

 26: /* Gauss point based evaluation */
 27: typedef struct {
 28:   PetscScalar gp_coords[NSD*GAUSS_POINTS];
 29:   PetscScalar eta[GAUSS_POINTS];
 30:   PetscScalar fx[GAUSS_POINTS];
 31:   PetscScalar fy[GAUSS_POINTS];
 32:   PetscScalar fz[GAUSS_POINTS];
 33:   PetscScalar hc[GAUSS_POINTS];
 34: } GaussPointCoefficients;

 36: typedef struct {
 37:   PetscScalar u_dof;
 38:   PetscScalar v_dof;
 39:   PetscScalar w_dof;
 40:   PetscScalar p_dof;
 41: } StokesDOF;

 43: typedef struct _p_CellProperties *CellProperties;
 44: struct _p_CellProperties {
 45:   PetscInt ncells;
 46:   PetscInt mx,my,mz;
 47:   PetscInt sex,sey,sez;
 48:   GaussPointCoefficients *gpc;
 49: };

 51: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz);

 53: /* elements */
 56: PetscErrorCode CellPropertiesCreate(DM da_stokes,CellProperties *C)
 57: {
 59:   CellProperties cells;
 60:   PetscInt mx,my,mz,sex,sey,sez;
 61: 
 63:   PetscMalloc(sizeof(struct _p_CellProperties),&cells);
 64: 
 65:   DMDAGetElementCorners(da_stokes,&sex,&sey,&sez,&mx,&my,&mz);
 66:   cells->mx = mx;
 67:   cells->my = my;
 68:   cells->mz = mz;
 69:   cells->ncells = mx * my * mz;
 70:   cells->sex = sex;
 71:   cells->sey = sey;
 72:   cells->sez = sez;
 73:   PetscMalloc(sizeof(GaussPointCoefficients)*mx*my*mz,&cells->gpc);
 74: 
 75:   *C = cells;
 76:   return(0);
 77: }

 81: PetscErrorCode CellPropertiesDestroy(CellProperties *C)
 82: {
 84:   CellProperties cells;

 87:   if (!C) { return(0); }
 88:   cells = *C;
 89:   PetscFree(cells->gpc);
 90:   PetscFree(cells);
 91:   *C = PETSC_NULL;
 92:   return(0);
 93: }

 97: PetscErrorCode CellPropertiesGetCell(CellProperties C,PetscInt I,PetscInt J,PetscInt K,GaussPointCoefficients **G)
 98: {
100:   *G = &C->gpc[ (I-C->sex) + (J-C->sey)*C->mx + (K-C->sez)*C->mx*C->my ];
101:   return(0);
102: }

104: /* FEM routines */
105: /*
106:  Element: Local basis function ordering
107:  1-----2
108:  |     |
109:  |     |
110:  0-----3
111:  */
112: static void ShapeFunctionQ13D_Evaluate(PetscScalar _xi[],PetscScalar Ni[])
113: {
114:   PetscReal xi   = PetscRealPart(_xi[0]);
115:   PetscReal eta  = PetscRealPart(_xi[1]);
116:   PetscReal zeta = PetscRealPart(_xi[2]);

118:   Ni[0] = 0.125 * ( 1.0 - xi ) * ( 1.0 - eta ) * ( 1.0 - zeta );
119:   Ni[1] = 0.125 * ( 1.0 - xi ) * ( 1.0 + eta ) * ( 1.0 - zeta );
120:   Ni[2] = 0.125 * ( 1.0 + xi ) * ( 1.0 + eta ) * ( 1.0 - zeta );
121:   Ni[3] = 0.125 * ( 1.0 + xi ) * ( 1.0 - eta ) * ( 1.0 - zeta );

123:   Ni[4] = 0.125 * ( 1.0 - xi ) * ( 1.0 - eta ) * ( 1.0 + zeta );
124:   Ni[5] = 0.125 * ( 1.0 - xi ) * ( 1.0 + eta ) * ( 1.0 + zeta );
125:   Ni[6] = 0.125 * ( 1.0 + xi ) * ( 1.0 + eta ) * ( 1.0 + zeta );
126:   Ni[7] = 0.125 * ( 1.0 + xi ) * ( 1.0 - eta ) * ( 1.0 + zeta );
127: }

129: static void ShapeFunctionQ13D_Evaluate_dxi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
130: {
131:   PetscReal xi   = PetscRealPart(_xi[0]);
132:   PetscReal eta  = PetscRealPart(_xi[1]);
133:   PetscReal zeta = PetscRealPart(_xi[2]);
134:   /* xi deriv */
135:   GNi[0][0] = - 0.125 * ( 1.0 - eta ) * ( 1.0 - zeta );
136:   GNi[0][1] = - 0.125 * ( 1.0 + eta ) * ( 1.0 - zeta );
137:   GNi[0][2] =   0.125 * ( 1.0 + eta ) * ( 1.0 - zeta );
138:   GNi[0][3] =   0.125 * ( 1.0 - eta ) * ( 1.0 - zeta );

140:   GNi[0][4] = - 0.125 * ( 1.0 - eta ) * ( 1.0 + zeta );
141:   GNi[0][5] = - 0.125 * ( 1.0 + eta ) * ( 1.0 + zeta );
142:   GNi[0][6] =   0.125 * ( 1.0 + eta ) * ( 1.0 + zeta );
143:   GNi[0][7] =   0.125 * ( 1.0 - eta ) * ( 1.0 + zeta );
144:   /* eta deriv */
145:   GNi[1][0] = - 0.125 * ( 1.0 - xi ) * ( 1.0 - zeta );
146:   GNi[1][1] =   0.125 * ( 1.0 - xi ) * ( 1.0 - zeta );
147:   GNi[1][2] =   0.125 * ( 1.0 + xi ) * ( 1.0 - zeta );
148:   GNi[1][3] = - 0.125 * ( 1.0 + xi ) * ( 1.0 - zeta );

150:   GNi[1][4] = - 0.125 * ( 1.0 - xi ) * ( 1.0 + zeta );
151:   GNi[1][5] =   0.125 * ( 1.0 - xi ) * ( 1.0 + zeta );
152:   GNi[1][6] =   0.125 * ( 1.0 + xi ) * ( 1.0 + zeta );
153:   GNi[1][7] = - 0.125 * ( 1.0 + xi ) * ( 1.0 + zeta );
154:   /* zeta deriv */
155:   GNi[2][0] = -0.125 * ( 1.0 - xi ) * ( 1.0 - eta );
156:   GNi[2][1] = -0.125 * ( 1.0 - xi ) * ( 1.0 + eta );
157:   GNi[2][2] = -0.125 * ( 1.0 + xi ) * ( 1.0 + eta );
158:   GNi[2][3] = -0.125 * ( 1.0 + xi ) * ( 1.0 - eta );

160:   GNi[2][4] = 0.125 * ( 1.0 - xi ) * ( 1.0 - eta );
161:   GNi[2][5] = 0.125 * ( 1.0 - xi ) * ( 1.0 + eta );
162:   GNi[2][6] = 0.125 * ( 1.0 + xi ) * ( 1.0 + eta );
163:   GNi[2][7] = 0.125 * ( 1.0 + xi ) * ( 1.0 - eta );
164: }

166: static void matrix_inverse_3x3(PetscScalar A[3][3],PetscScalar B[3][3])
167: {
168:   PetscScalar t4, t6, t8, t10, t12, t14, t17;

170:   t4 = A[2][0] * A[0][1];
171:   t6 = A[2][0] * A[0][2];
172:   t8 = A[1][0] * A[0][1];
173:   t10 = A[1][0] * A[0][2];
174:   t12 = A[0][0] * A[1][1];
175:   t14 = A[0][0] * A[1][2];
176:   t17 = 0.1e1 / (t4 * A[1][2] - t6 * A[1][1] - t8 * A[2][2] + t10 * A[2][1] + t12 * A[2][2] - t14 * A[2][1]);

178:   B[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) * t17;
179:   B[0][1] = -(A[0][1] * A[2][2] - A[0][2] * A[2][1]) * t17;
180:   B[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) * t17;
181:   B[1][0] = -(-A[2][0] * A[1][2] + A[1][0] * A[2][2]) * t17;
182:   B[1][1] = (-t6 + A[0][0] * A[2][2]) * t17;
183:   B[1][2] = -(-t10 + t14) * t17;
184:   B[2][0] = (-A[2][0] * A[1][1] + A[1][0] * A[2][1]) * t17;
185:   B[2][1] = -(-t4 + A[0][0] * A[2][1]) * t17;
186:   B[2][2] = (-t8 + t12) * t17;
187: }

189: static void ShapeFunctionQ13D_Evaluate_dx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
190: {
191:   PetscScalar J00,J01,J02,J10,J11,J12,J20,J21,J22;
192:   PetscInt n;
193:   PetscScalar iJ[3][3],JJ[3][3];

195:   J00 = J01 = J02 = 0.0;
196:   J10 = J11 = J12 = 0.0;
197:   J20 = J21 = J22 = 0.0;
198:   for (n=0; n<NODES_PER_EL; n++) {
199:     PetscScalar cx = coords[ NSD*n + 0 ];
200:     PetscScalar cy = coords[ NSD*n + 1 ];
201:     PetscScalar cz = coords[ NSD*n + 2 ];

203:     /* J_ij = d(x_j) / d(xi_i) *//* J_ij = \sum _I GNi[j][I} * x_i */
204:     J00 = J00 + GNi[0][n] * cx;  /* J_xx */
205:     J01 = J01 + GNi[0][n] * cy;  /* J_xy = dx/deta */
206:     J02 = J02 + GNi[0][n] * cz;  /* J_xz = dx/dzeta */

208:     J10 = J10 + GNi[1][n] * cx;  /* J_yx = dy/dxi */
209:     J11 = J11 + GNi[1][n] * cy;  /* J_yy */
210:     J12 = J12 + GNi[1][n] * cz;  /* J_yz */

212:     J20 = J20 + GNi[2][n] * cx;  /* J_zx */
213:     J21 = J21 + GNi[2][n] * cy;  /* J_zy */
214:     J22 = J22 + GNi[2][n] * cz;  /* J_zz */
215:   }

217:   JJ[0][0] = J00;                        JJ[0][1] = J01;                        JJ[0][2] = J02;
218:   JJ[1][0] = J10;                        JJ[1][1] = J11;                        JJ[1][2] = J12;
219:   JJ[2][0] = J20;                        JJ[2][1] = J21;                        JJ[2][2] = J22;

221:   matrix_inverse_3x3(JJ,iJ);

223:   *det_J = J00*J11*J22 - J00*J12*J21 - J10*J01*J22 + J10*J02*J21 + J20*J01*J12 - J20*J02*J11;

225:   for (n=0; n<NODES_PER_EL; n++) {
226:     GNx[0][n] = GNi[0][n]*iJ[0][0] + GNi[1][n]*iJ[0][1] + GNi[2][n]*iJ[0][2];
227:     GNx[1][n] = GNi[0][n]*iJ[1][0] + GNi[1][n]*iJ[1][1] + GNi[2][n]*iJ[1][2];
228:     GNx[2][n] = GNi[0][n]*iJ[2][0] + GNi[1][n]*iJ[2][1] + GNi[2][n]*iJ[2][2];
229:   }
230: }

232: static void ConstructGaussQuadrature3D(PetscInt *ngp,PetscScalar gp_xi[][NSD],PetscScalar gp_weight[])
233: {
234:   *ngp         = 8;
235:   gp_xi[0][0]  = -0.57735026919; gp_xi[0][1] = -0.57735026919; gp_xi[0][2] = -0.57735026919;
236:   gp_xi[1][0]  = -0.57735026919; gp_xi[1][1] =  0.57735026919; gp_xi[1][2] = -0.57735026919;
237:   gp_xi[2][0]  =  0.57735026919; gp_xi[2][1] =  0.57735026919; gp_xi[2][2] = -0.57735026919;
238:   gp_xi[3][0]  =  0.57735026919; gp_xi[3][1] = -0.57735026919; gp_xi[3][2] = -0.57735026919;

240:   gp_xi[4][0]  = -0.57735026919; gp_xi[4][1] = -0.57735026919; gp_xi[4][2] =  0.57735026919;
241:   gp_xi[5][0]  = -0.57735026919; gp_xi[5][1] =  0.57735026919; gp_xi[5][2] =  0.57735026919;
242:   gp_xi[6][0]  =  0.57735026919; gp_xi[6][1] =  0.57735026919; gp_xi[6][2] =  0.57735026919;
243:   gp_xi[7][0]  =  0.57735026919; gp_xi[7][1] = -0.57735026919; gp_xi[7][2] =  0.57735026919;

245:   gp_weight[0] = 1.0;
246:   gp_weight[1] = 1.0;
247:   gp_weight[2] = 1.0;
248:   gp_weight[3] = 1.0;

250:   gp_weight[4] = 1.0;
251:   gp_weight[5] = 1.0;
252:   gp_weight[6] = 1.0;
253:   gp_weight[7] = 1.0;
254: }

256: /* procs to the left claim the ghost node as their element */
259: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
260: {
261:   PetscInt m,n,p,M,N,P;
262:   PetscInt sx,sy,sz;

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

269:   if (mxl) {
270:     *mxl = m;
271:     if ((sx+m) == M) {  /* last proc */
272:       *mxl = m-1;
273:     }
274:   }
275:   if (myl) {
276:     *myl = n;
277:     if ((sy+n) == N) {  /* last proc */
278:       *myl = n-1;
279:     }
280:   }
281:   if (mzl) {
282:     *mzl = p;
283:     if ((sz+p) == P) {  /* last proc */
284:       *mzl = p-1;
285:     }
286:   }
287:   return(0);
288: }

292: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
293: {
294:   PetscInt si,sj,sk;

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

300:   if (sx) {
301:     *sx = si;
302:     if (si != 0) {
303:       *sx = si+1;
304:     }
305:   }
306:   if (sy) {
307:     *sy = sj;
308:     if (sj != 0) {
309:       *sy = sj+1;
310:     }
311:   }
312:   if (sz) {
313:     *sz = sk;
314:     if (sk != 0) {
315:       *sz = sk+1;
316:     }
317:   }
318:   DMDAGetLocalElementSize(da,mx,my,mz);
319:   return(0);
320: }

322: /*
323:  i,j are the element indices
324:  The unknown is a vector quantity.
325:  The s[].c is used to indicate the degree of freedom.
326:  */
329: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j,PetscInt k)
330: {
331:   PetscInt n;
333:   /* velocity */
334:   n = 0;
335:   /* node 0 */
336:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++; /* Vx0 */
337:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++; /* Vy0 */
338:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++; /* Vz0 */

340:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
341:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
342:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;

344:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
345:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
346:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;

348:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++;
349:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++;
350:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++;

352:   /* */
353:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++; /* Vx4 */
354:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++; /* Vy4 */
355:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++; /* Vz4 */

357:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
358:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
359:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;

361:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
362:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
363:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;

365:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++;
366:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++;
367:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++;

369:   /* pressure */
370:   n = 0;
371:   s_p[n].i = i;   s_p[n].j = j;   s_p[n].k = k; s_p[n].c = 3; n++; /* P0 */
372:   s_p[n].i = i;   s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
373:   s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
374:   s_p[n].i = i+1; s_p[n].j = j;   s_p[n].k = k; s_p[n].c = 3; n++;

376:   s_p[n].i = i;   s_p[n].j = j;   s_p[n].k = k+1; s_p[n].c = 3; n++; /* P0 */
377:   s_p[n].i = i;   s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
378:   s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
379:   s_p[n].i = i+1; s_p[n].j = j;   s_p[n].k = k+1; s_p[n].c = 3; n++;
380:   return(0);
381: }

385: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords,PetscInt i,PetscInt j,PetscInt k,PetscScalar el_coord[])
386: {
388:   /* get coords for the element */
389:   el_coord[0] = coords[k  ][j  ][i  ].x;
390:   el_coord[1] = coords[k  ][j  ][i  ].y;
391:   el_coord[2] = coords[k  ][j  ][i  ].z;

393:   el_coord[3] = coords[k  ][j+1][i].x;
394:   el_coord[4] = coords[k  ][j+1][i].y;
395:   el_coord[5] = coords[k  ][j+1][i].z;

397:   el_coord[6] = coords[k  ][j+1][i+1].x;
398:   el_coord[7] = coords[k  ][j+1][i+1].y;
399:   el_coord[8] = coords[k  ][j+1][i+1].z;

401:   el_coord[9 ] = coords[k  ][j  ][i+1].x;
402:   el_coord[10] = coords[k  ][j  ][i+1].y;
403:   el_coord[11] = coords[k  ][j  ][i+1].z;

405:   el_coord[12] = coords[k+1][j  ][i  ].x;
406:   el_coord[13] = coords[k+1][j  ][i  ].y;
407:   el_coord[14] = coords[k+1][j  ][i  ].z;

409:   el_coord[15] = coords[k+1][j+1][i  ].x;
410:   el_coord[16] = coords[k+1][j+1][i  ].y;
411:   el_coord[17] = coords[k+1][j+1][i  ].z;

413:   el_coord[18] = coords[k+1][j+1][i+1].x;
414:   el_coord[19] = coords[k+1][j+1][i+1].y;
415:   el_coord[20] = coords[k+1][j+1][i+1].z;

417:   el_coord[21] = coords[k+1][j  ][i+1].x;
418:   el_coord[22] = coords[k+1][j  ][i+1].y;
419:   el_coord[23] = coords[k+1][j  ][i+1].z;
420:   return(0);
421: }

425: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field,PetscInt i,PetscInt j,PetscInt k,StokesDOF nodal_fields[])
426: {
428:   /* get the nodal fields for u */
429:   nodal_fields[0].u_dof = field[k  ][j  ][i  ].u_dof;
430:   nodal_fields[0].v_dof = field[k  ][j  ][i  ].v_dof;
431:   nodal_fields[0].w_dof = field[k  ][j  ][i  ].w_dof;

433:   nodal_fields[1].u_dof = field[k  ][j+1][i  ].u_dof;
434:   nodal_fields[1].v_dof = field[k  ][j+1][i  ].v_dof;
435:   nodal_fields[1].w_dof = field[k  ][j+1][i  ].w_dof;

437:   nodal_fields[2].u_dof = field[k  ][j+1][i+1].u_dof;
438:   nodal_fields[2].v_dof = field[k  ][j+1][i+1].v_dof;
439:   nodal_fields[2].w_dof = field[k  ][j+1][i+1].w_dof;

441:   nodal_fields[3].u_dof = field[k  ][j  ][i+1].u_dof;
442:   nodal_fields[3].v_dof = field[k  ][j  ][i+1].v_dof;
443:   nodal_fields[3].w_dof = field[k  ][j  ][i+1].w_dof;

445:   nodal_fields[4].u_dof = field[k+1][j  ][i  ].u_dof;
446:   nodal_fields[4].v_dof = field[k+1][j  ][i  ].v_dof;
447:   nodal_fields[4].w_dof = field[k+1][j  ][i  ].w_dof;

449:   nodal_fields[5].u_dof = field[k+1][j+1][i  ].u_dof;
450:   nodal_fields[5].v_dof = field[k+1][j+1][i  ].v_dof;
451:   nodal_fields[5].w_dof = field[k+1][j+1][i  ].w_dof;

453:   nodal_fields[6].u_dof = field[k+1][j+1][i+1].u_dof;
454:   nodal_fields[6].v_dof = field[k+1][j+1][i+1].v_dof;
455:   nodal_fields[6].w_dof = field[k+1][j+1][i+1].w_dof;

457:   nodal_fields[7].u_dof = field[k+1][j  ][i+1].u_dof;
458:   nodal_fields[7].v_dof = field[k+1][j  ][i+1].v_dof;
459:   nodal_fields[7].w_dof = field[k+1][j  ][i+1].w_dof;

461:   /* get the nodal fields for p */
462:   nodal_fields[0].p_dof = field[k  ][j  ][i  ].p_dof;
463:   nodal_fields[1].p_dof = field[k  ][j+1][i  ].p_dof;
464:   nodal_fields[2].p_dof = field[k  ][j+1][i+1].p_dof;
465:   nodal_fields[3].p_dof = field[k  ][j  ][i+1].p_dof;

467:   nodal_fields[4].p_dof = field[k+1][j  ][i  ].p_dof;
468:   nodal_fields[5].p_dof = field[k+1][j+1][i  ].p_dof;
469:   nodal_fields[6].p_dof = field[k+1][j+1][i+1].p_dof;
470:   nodal_fields[7].p_dof = field[k+1][j  ][i+1].p_dof;
471:   return(0);
472: }

474: static PetscInt ASS_MAP_wIwDI_uJuDJ(
475:   PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,
476:   PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)
477: {
478:   PetscInt ij;
479:   PETSC_UNUSED PetscInt r,c,nr,nc;

481:   nr = w_NPE*w_dof;
482:   nc = u_NPE*u_dof;

484:   r = w_dof*wi+wd;
485:   c = u_dof*ui+ud;

487:   ij = r*nc+c;

489:   return ij;
490: }

494: static PetscErrorCode DMDASetValuesLocalStencil3D_ADD_VALUES(StokesDOF ***fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
495: {
496:   PetscInt n,I,J,K;

499:   for (n = 0; n<NODES_PER_EL; n++) {
500:     I = u_eqn[NSD*n  ].i;
501:     J = u_eqn[NSD*n  ].j;
502:     K = u_eqn[NSD*n  ].k;
503:     fields_F[K][J][I].u_dof = fields_F[K][J][I].u_dof+Fe_u[NSD*n  ];

505:     I = u_eqn[NSD*n+1].i;
506:     J = u_eqn[NSD*n+1].j;
507:     K = u_eqn[NSD*n+1].k;
508:     fields_F[K][J][I].v_dof = fields_F[K][J][I].v_dof+Fe_u[NSD*n+1];

510:     I = u_eqn[NSD*n+2].i;
511:     J = u_eqn[NSD*n+2].j;
512:     K = u_eqn[NSD*n+2].k;
513:     fields_F[K][J][I].w_dof = fields_F[K][J][I].w_dof+Fe_u[NSD*n+2];

515:     I = p_eqn[n].i;
516:     J = p_eqn[n].j;
517:     K = p_eqn[n].k;
518:     fields_F[K][J][I].p_dof = fields_F[K][J][I].p_dof+Fe_p[n ];

520:   }
521:   return(0);
522: }

524: static void FormStressOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
525: {
526:   PetscInt       ngp;
527:   PetscScalar    gp_xi[GAUSS_POINTS][NSD];
528:   PetscScalar    gp_weight[GAUSS_POINTS];
529:   PetscInt       p,i,j,k;
530:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
531:   PetscScalar    J_p,tildeD[6];
532:   PetscScalar    B[6][U_DOFS*NODES_PER_EL];
533:   const PetscInt nvdof = U_DOFS*NODES_PER_EL;

535:   /* define quadrature rule */
536:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

538:   /* evaluate integral */
539:   for (p = 0; p < ngp; p++) {
540:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
541:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);

543:     for (i = 0; i < NODES_PER_EL; i++) {
544:       PetscScalar d_dx_i = GNx_p[0][i];
545:       PetscScalar d_dy_i = GNx_p[1][i];
546:       PetscScalar d_dz_i = GNx_p[2][i];

548:       B[0][3*i  ] = d_dx_i; B[0][3*i+1] = 0.0;    B[0][3*i+2] = 0.0;
549:       B[1][3*i  ] = 0.0;    B[1][3*i+1] = d_dy_i; B[1][3*i+2] = 0.0;
550:       B[2][3*i  ] = 0.0;    B[2][3*i+1] = 0.0;    B[2][3*i+2] = d_dz_i;

552:       B[3][3*i] = d_dy_i; B[3][3*i+1] = d_dx_i;  B[3][3*i+2] = 0.0;   /* e_xy */
553:       B[4][3*i] = d_dz_i; B[4][3*i+1] = 0.0;     B[4][3*i+2] = d_dx_i;/* e_xz */
554:       B[5][3*i] = 0.0;    B[5][3*i+1] = d_dz_i;  B[5][3*i+2] = d_dy_i;/* e_yz */
555:     }


558:     tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
559:     tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
560:     tildeD[2] = 2.0*gp_weight[p]*J_p*eta[p];

562:     tildeD[3] =     gp_weight[p]*J_p*eta[p];
563:     tildeD[4] =     gp_weight[p]*J_p*eta[p];
564:     tildeD[5] =     gp_weight[p]*J_p*eta[p];

566:     /* form Bt tildeD B */
567:     /*
568:      Ke_ij = Bt_ik . D_kl . B_lj
569:      = B_ki . D_kl . B_lj
570:      = B_ki . D_kk . B_kj
571:      */
572:     for (i = 0; i < nvdof; i++) {
573:       for (j = i; j < nvdof; j++) {
574:         for (k = 0; k < 6; k++) {
575:           Ke[i*nvdof+j] += B[k][i]*tildeD[k]*B[k][j];
576:         }
577:       }
578:     }

580:   }
581:   /* fill lower triangular part */
582: #ifdef ASSEMBLE_LOWER_TRIANGULAR
583:   for (i = 0; i < nvdof; i++) {
584:     for (j = i; j < nvdof; j++) {
585:       Ke[j*nvdof+i] = Ke[i*nvdof+j];
586:     }
587:   }
588: #endif
589: }

591: static void FormGradientOperatorQ13D(PetscScalar Ke[],PetscScalar coords[])
592: {
593:   PetscInt    ngp;
594:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
595:   PetscScalar gp_weight[GAUSS_POINTS];
596:   PetscInt    p,i,j,di;
597:   PetscScalar Ni_p[NODES_PER_EL];
598:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
599:   PetscScalar J_p,fac;

601:   /* define quadrature rule */
602:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

604:   /* evaluate integral */
605:   for (p = 0; p < ngp; p++) {
606:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
607:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
608:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
609:     fac = gp_weight[p]*J_p;

611:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
612:       for (di = 0; di < NSD; di++) { /* u dofs */
613:         for (j = 0; j < NODES_PER_EL; j++) {  /* p nodes, p dofs = 1 (ie no loop) */
614:           PetscInt IJ;
615:           IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,3,j,0,NODES_PER_EL,1);

617:           Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
618:         }
619:       }
620:     }
621:   }
622: }

624: static void FormDivergenceOperatorQ13D(PetscScalar De[],PetscScalar coords[])
625: {
626:   PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
627:   PetscInt    i,j;
628:   PetscInt    nr_g,nc_g;

630:   PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
631:   FormGradientOperatorQ13D(Ge,coords);

633:   nr_g = U_DOFS*NODES_PER_EL;
634:   nc_g = P_DOFS*NODES_PER_EL;

636:   for (i = 0; i < nr_g; i++) {
637:     for (j = 0; j < nc_g; j++) {
638:       De[nr_g*j+i] = Ge[nc_g*i+j];
639:     }
640:   }
641: }

643: static void FormStabilisationOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
644: {
645:   PetscInt    ngp;
646:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
647:   PetscScalar gp_weight[GAUSS_POINTS];
648:   PetscInt    p,i,j;
649:   PetscScalar Ni_p[NODES_PER_EL];
650:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
651:   PetscScalar J_p,fac,eta_avg;

653:   /* define quadrature rule */
654:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

656:   /* evaluate integral */
657:   for (p = 0; p < ngp; p++) {
658:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
659:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
660:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
661:     fac = gp_weight[p]*J_p;
662:     /*
663:      for (i = 0; i < NODES_PER_EL; i++) {
664:      for (j = i; j < NODES_PER_EL; j++) {
665:      Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
666:      }
667:      }
668:      */

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

677:   /* scale */
678:   eta_avg = 0.0;
679:   for (p = 0; p < ngp; p++) {
680:     eta_avg += eta[p];
681:   }
682:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
683:   fac     = 1.0/eta_avg;

685:   /*
686:    for (i = 0; i < NODES_PER_EL; i++) {
687:    for (j = i; j < NODES_PER_EL; j++) {
688:    Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
689:    #ifdef ASSEMBLE_LOWER_TRIANGULAR
690:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
691:    #endif
692:    }
693:    }
694:    */

696:   for (i = 0; i < NODES_PER_EL; i++) {
697:     for (j = 0; j < NODES_PER_EL; j++) {
698:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
699:     }
700:   }
701: }

703: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
704: {
705:   PetscInt    ngp;
706:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
707:   PetscScalar gp_weight[GAUSS_POINTS];
708:   PetscInt    p,i,j;
709:   PetscScalar Ni_p[NODES_PER_EL];
710:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
711:   PetscScalar J_p,fac,eta_avg;

713:   /* define quadrature rule */
714:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

716:   /* evaluate integral */
717:   for (p = 0; p < ngp; p++) {
718:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
719:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
720:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
721:     fac = gp_weight[p]*J_p;

723:     /*
724:      for (i = 0; i < NODES_PER_EL; i++) {
725:      for (j = i; j < NODES_PER_EL; j++) {
726:      Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
727:      }
728:      }
729:      */

731:     for (i = 0; i < NODES_PER_EL; i++) {
732:       for (j = 0; j < NODES_PER_EL; j++) {
733:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
734:       }
735:     }
736:   }

738:   /* scale */
739:   eta_avg = 0.0;
740:   for (p = 0; p < ngp; p++) {
741:     eta_avg += eta[p];
742:   }
743:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
744:   fac     = 1.0/eta_avg;
745:   /*
746:    for (i = 0; i < NODES_PER_EL; i++) {
747:    for (j = i; j < NODES_PER_EL; j++) {
748:    Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
749:    #ifdef ASSEMBLE_LOWER_TRIANGULAR
750:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
751:    #endif
752:    }
753:    }
754:    */

756:   for (i = 0; i < NODES_PER_EL; i++) {
757:     for (j = 0; j < NODES_PER_EL; j++) {
758:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
759:     }
760:   }
761: }

763: static void FormMomentumRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[],PetscScalar fz[])
764: {
765:   PetscInt    ngp;
766:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
767:   PetscScalar gp_weight[GAUSS_POINTS];
768:   PetscInt    p,i;
769:   PetscScalar Ni_p[NODES_PER_EL];
770:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
771:   PetscScalar J_p,fac;

773:   /* define quadrature rule */
774:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

776:   /* evaluate integral */
777:   for (p = 0; p < ngp; p++) {
778:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
779:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
780:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
781:     fac = gp_weight[p]*J_p;

783:     for (i = 0; i < NODES_PER_EL; i++) {
784:       Fe[NSD*i  ] -= fac*Ni_p[i]*fx[p];
785:       Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
786:       Fe[NSD*i+2] -= fac*Ni_p[i]*fz[p];
787:     }
788:   }
789: }

791: static void FormContinuityRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar hc[])
792: {
793:   PetscInt    ngp;
794:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
795:   PetscScalar gp_weight[GAUSS_POINTS];
796:   PetscInt    p,i;
797:   PetscScalar Ni_p[NODES_PER_EL];
798:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
799:   PetscScalar J_p,fac;

801:   /* define quadrature rule */
802:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

804:   /* evaluate integral */
805:   for (p = 0; p < ngp; p++) {
806:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
807:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
808:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
809:     fac = gp_weight[p]*J_p;

811:     for (i = 0; i < NODES_PER_EL; i++) {
812:       Fe[i] -= fac*Ni_p[i]*hc[p];
813:     }
814:   }
815: }

817: #define _ZERO_ROWCOL_i(A,i) {                   \
818:     PetscInt KK;                                \
819:     PetscScalar tmp = A[24*(i)+(i)];            \
820:     for(KK=0;KK<24;KK++){A[24*(i)+KK]=0.0;}     \
821:     for(KK=0;KK<24;KK++){A[24*KK+(i)]=0.0;}     \
822:     A[24*(i)+(i)] = tmp;}                       \

824: #define _ZERO_ROW_i(A,i) {                      \
825:     PetscInt KK;                                \
826:     for(KK=0;KK<8;KK++){A[8*(i)+KK]=0.0;}}

828: #define _ZERO_COL_i(A,i) {                      \
829:     PetscInt KK;                                \
830:     for(KK=0;KK<8;KK++){A[24*KK+(i)]=0.0;}}

834: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,CellProperties cell_properties)
835: {
836:   DM                     cda;
837:   Vec                    coords;
838:   DMDACoor3d             ***_coords;
839:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS];
840:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS];
841:   PetscInt               sex,sey,sez,mx,my,mz;
842:   PetscInt               ei,ej,ek;
843:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
844:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
845:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
846:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
847:   PetscScalar            el_coords[NODES_PER_EL*NSD];
848:   GaussPointCoefficients *props;
849:   PetscScalar            *prop_eta;
850:   PetscInt               n,M,N,P;
851:   PetscLogDouble         t0,t1;
852:   PetscErrorCode         ierr;


856:   DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
857:   /* setup for coords */
858:   DMDAGetCoordinateDA(stokes_da,&cda);
859:   DMDAGetGhostedCoordinates(stokes_da,&coords);
860:   DMDAVecGetArray(cda,coords,&_coords);

862:   DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
863:   PetscGetTime(&t0);
864:   for (ek = sez; ek < sez+mz; ek++) {
865:     for (ej = sey; ej < sey+my; ej++) {
866:       for (ei = sex; ei < sex+mx; ei++) {
867:         /* get coords for the element */
868:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
869:         /* get cell properties */
870:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
871:         /* get coefficients for the element */
872:         prop_eta = props->eta;

874:         /* initialise element stiffness matrix */
875:         PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
876:         PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
877:         PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
878:         PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);

880:         /* form element stiffness matrix */
881:         FormStressOperatorQ13D(Ae,el_coords,prop_eta);
882:         FormGradientOperatorQ13D(Ge,el_coords);
883:         /*#ifdef ASSEMBLE_LOWER_TRIANGULAR*/
884:         FormDivergenceOperatorQ13D(De,el_coords);
885:         /*#endif*/
886:         FormStabilisationOperatorQ13D(Ce,el_coords,prop_eta);

888:         /* insert element matrix into global matrix */
889:         DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);

891:         for (n=0; n<NODES_PER_EL; n++) {
892:           if ( (u_eqn[3*n  ].i == 0) || (u_eqn[3*n  ].i == M-1) ) {
893:             _ZERO_ROWCOL_i(Ae,3*n);
894:             _ZERO_ROW_i   (Ge,3*n);
895:             _ZERO_COL_i   (De,3*n);
896:           }

898:           if ( (u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1) ) {
899:             _ZERO_ROWCOL_i(Ae,3*n+1);
900:             _ZERO_ROW_i   (Ge,3*n+1);
901:             _ZERO_COL_i   (De,3*n+1);
902:           }

904:           if ( (u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1) ) {
905:             _ZERO_ROWCOL_i(Ae,3*n+2);
906:             _ZERO_ROW_i   (Ge,3*n+2);
907:             _ZERO_COL_i   (De,3*n+2);
908:           }
909:         }
910:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
911:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
912:         MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
913:         MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
914:       }
915:     }
916:   }
917:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
918:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

922:   PetscGetTime(&t1);
923:   return(0);
924: }

928: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,CellProperties cell_properties)
929: {
930:   DM                     cda;
931:   Vec                    coords;
932:   DMDACoor3d             ***_coords;
933:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS];
934:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS];
935:   PetscInt               sex,sey,sez,mx,my,mz;
936:   PetscInt               ei,ej,ek;
937:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
938:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
939:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
940:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
941:   PetscScalar            el_coords[NODES_PER_EL*NSD];
942:   GaussPointCoefficients *props;
943:   PetscScalar            *prop_eta;
944:   PetscInt               n,M,N,P;
945:   PetscErrorCode         ierr;


949:   DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
950:   /* setup for coords */
951:   DMDAGetCoordinateDA(stokes_da,&cda);
952:   DMDAGetGhostedCoordinates(stokes_da,&coords);
953:   DMDAVecGetArray(cda,coords,&_coords);

955:   DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
956:   for (ek = sez; ek < sez+mz; ek++) {
957:     for (ej = sey; ej < sey+my; ej++) {
958:       for (ei = sex; ei < sex+mx; ei++) {
959:         /* get coords for the element */
960:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
961:         /* get cell properties */
962:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
963:         /* get coefficients for the element */
964:         prop_eta = props->eta;

966:         /* initialise element stiffness matrix */
967:         PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
968:         PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
969:         PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
970:         PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);

972:         /* form element stiffness matrix */
973:         FormStressOperatorQ13D(Ae,el_coords,prop_eta);
974:         FormGradientOperatorQ13D(Ge,el_coords);
975:         /* FormDivergenceOperatorQ13D(De,el_coords); */
976:         FormScaledMassMatrixOperatorQ13D(Ce,el_coords,prop_eta);

978:         /* insert element matrix into global matrix */
979:         DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);

981:         for (n=0; n<NODES_PER_EL; n++) {
982:           if ( (u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1) ) {
983:             _ZERO_ROWCOL_i(Ae,3*n);
984:             _ZERO_ROW_i   (Ge,3*n);
985:             _ZERO_COL_i   (De,3*n);
986:           }

988:           if ( (u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1) ) {
989:             _ZERO_ROWCOL_i(Ae,3*n+1);
990:             _ZERO_ROW_i   (Ge,3*n+1);
991:             _ZERO_COL_i   (De,3*n+1);
992:           }

994:           if ( (u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1) ) {
995:             _ZERO_ROWCOL_i(Ae,3*n+2);
996:             _ZERO_ROW_i   (Ge,3*n+2);
997:             _ZERO_COL_i   (De,3*n+2);
998:           }
999:         }
1000:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
1001:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
1002:         /*MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);*/
1003:         MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
1004:       }
1005:     }
1006:   }
1007:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1008:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

1010:   DMDAVecRestoreArray(cda,coords,&_coords);
1011:   return(0);
1012: }

1016: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,CellProperties cell_properties)
1017: {
1018:   DM                     cda;
1019:   Vec                    coords;
1020:   DMDACoor3d             ***_coords;
1021:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS];
1022:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS];
1023:   PetscInt               sex,sey,sez,mx,my,mz;
1024:   PetscInt               ei,ej,ek;
1025:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
1026:   PetscScalar            He[NODES_PER_EL*P_DOFS];
1027:   PetscScalar            el_coords[NODES_PER_EL*NSD];
1028:   GaussPointCoefficients *props;
1029:   PetscScalar            *prop_fx,*prop_fy,*prop_fz,*prop_hc;
1030:   Vec                    local_F;
1031:   StokesDOF              ***ff;
1032:   PetscInt               n,M,N,P;
1033:   PetscErrorCode         ierr;


1037:   DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
1038:   /* setup for coords */
1039:   DMDAGetCoordinateDA(stokes_da,&cda);
1040:   DMDAGetGhostedCoordinates(stokes_da,&coords);
1041:   DMDAVecGetArray(cda,coords,&_coords);

1043:   /* get acces to the vector */
1044:   DMGetLocalVector(stokes_da,&local_F);
1045:   VecZeroEntries(local_F);
1046:   DMDAVecGetArray(stokes_da,local_F,&ff);

1048:   DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
1049:   for (ek = sez; ek < sez+mz; ek++) {
1050:     for (ej = sey; ej < sey+my; ej++) {
1051:       for (ei = sex; ei < sex+mx; ei++) {
1052:         /* get coords for the element */
1053:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1054:         /* get cell properties */
1055:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
1056:         /* get coefficients for the element */
1057:         prop_fx = props->fx;
1058:         prop_fy = props->fy;
1059:         prop_fz = props->fz;
1060:         prop_hc = props->hc;

1062:         /* initialise element stiffness matrix */
1063:         PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
1064:         PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);

1066:         /* form element stiffness matrix */
1067:         FormMomentumRhsQ13D(Fe,el_coords,prop_fx,prop_fy,prop_fz);
1068:         FormContinuityRhsQ13D(He,el_coords,prop_hc);

1070:         /* insert element matrix into global matrix */
1071:         DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);

1073:         for (n=0; n<NODES_PER_EL; n++) {
1074:           if ( (u_eqn[3*n  ].i == 0) || (u_eqn[3*n  ].i == M-1) ) {
1075:             Fe[3*n  ] = 0.0;
1076:           }

1078:           if ( (u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1) ) {
1079:             Fe[3*n+1] = 0.0;
1080:           }

1082:           if ( (u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1) ) {
1083:             Fe[3*n+2] = 0.0;
1084:           }
1085:         }

1087:         DMDASetValuesLocalStencil3D_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
1088:       }
1089:     }
1090:   }
1091:   DMDAVecRestoreArray(stokes_da,local_F,&ff);
1092:   DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
1093:   DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
1094:   DMRestoreLocalVector(stokes_da,&local_F);

1096:   DMDAVecRestoreArray(cda,coords,&_coords);
1097:   return(0);
1098: }

1100: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta,PetscReal *MX,PetscReal *MY,PetscReal *MZ)
1101: {
1102:   *theta = 0.0;
1103:   *MX = 2.0 * M_PI;
1104:   *MY = 2.0 * M_PI;
1105:   *MZ = 2.0 * M_PI;
1106: }
1107: static void evaluate_MS_FrankKamentski(PetscReal pos[],PetscReal v[],PetscReal *p,PetscReal *eta,PetscReal Fm[],PetscReal *Fc)
1108: {
1109:   PetscReal x,y,z;
1110:   PetscReal theta,MX,MY,MZ;

1112:   evaluate_MS_FrankKamentski_constants(&theta,&MX,&MY,&MZ);
1113:   x = pos[0];
1114:   y = pos[1];
1115:   z = pos[2];
1116:   if (v) {
1117:     /*
1118:      v[0] = pow(z,3)*exp(y)*sin(M_PI*x);
1119:      v[1] = z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y);
1120:      v[2] = -(pow(x,3) + pow(y,3))*sin(2.0*M_PI*z);
1121:      */
1122:     /*
1123:      v[0] = sin(M_PI*x);
1124:      v[1] = sin(M_PI*y);
1125:      v[2] = sin(M_PI*z);
1126:      */
1127:     v[0] = pow(z,3)*exp(y)*sin(M_PI*x);
1128:     v[1] = z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y);
1129:     v[2] = pow(z,2)*(cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 - M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2) - M_PI*pow(z,4)*cos(M_PI*x)*exp(y)/4 ;
1130:   }
1131:   if (p) {
1132:     *p = pow(x,2) + pow(y,2) + pow(z,2);
1133:   }
1134:   if (eta) {
1135:     /**eta = exp(-theta*(1.0 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)));*/
1136:     *eta = 1.0;
1137:   }
1138:   if (Fm) {
1139:     /*
1140:      Fm[0] = -2*x - 2.0*pow(M_PI,2)*pow(z,3)*exp(y)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(M_PI*x) - 0.2*MZ*theta*(-1.5*pow(x,2)*sin(2.0*M_PI*z) + 1.5*pow(z,2)*exp(y)*sin(M_PI*x))*cos(MX*x)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MY*y)*sin(MZ*z) - 0.2*M_PI*MX*theta*pow(z,3)*cos(M_PI*x)*cos(MZ*z)*exp(y)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MX*x)*sin(MY*y) + 2.0*(3.0*z*exp(y)*sin(M_PI*x) - 3.0*M_PI*pow(x,2)*cos(2.0*M_PI*z))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*(0.5*pow(z,3)*exp(y)*sin(M_PI*x) + M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x) - 1.0*z*pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*theta*(1 + 0.1*MY*cos(MX*x)*cos(MY*y)*cos(MZ*z))*(0.5*pow(z,3)*exp(y)*sin(M_PI*x) - 1.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) ;
1141:      Fm[1] = -2*y - 0.2*MX*theta*(0.5*pow(z,3)*exp(y)*sin(M_PI*x) - 1.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x))*cos(MZ*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MX*x)*sin(MY*y) - 0.2*MZ*theta*(-1.5*pow(y,2)*sin(2.0*M_PI*z) + 0.5*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y))*cos(MX*x)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MY*y)*sin(MZ*z) + 2.0*(-2.0*z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 0.5*M_PI*pow(z,3)*cos(M_PI*x)*exp(y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*(z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 2*M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*theta*(1 + 0.1*MY*cos(MX*x)*cos(MY*y)*cos(MZ*z))*(-z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) - 6.0*M_PI*pow(y,2)*cos(2.0*M_PI*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)));
1142:      Fm[2] = -2*z + 8.0*pow(M_PI,2)*(pow(x,3) + pow(y,3))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(2.0*M_PI*z) - 0.2*MX*theta*(-1.5*pow(x,2)*sin(2.0*M_PI*z) + 1.5*pow(z,2)*exp(y)*sin(M_PI*x))*cos(MZ*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MX*x)*sin(MY*y) + 0.4*M_PI*MZ*theta*(pow(x,3) + pow(y,3))*cos(MX*x)*cos(2.0*M_PI*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MY*y)*sin(MZ*z) + 2.0*(-3.0*x*sin(2.0*M_PI*z) + 1.5*M_PI*pow(z,2)*cos(M_PI*x)*exp(y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*(-3.0*y*sin(2.0*M_PI*z) - 0.5*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 0.5*M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*theta*(1 + 0.1*MY*cos(MX*x)*cos(MY*y)*cos(MZ*z))*(-1.5*pow(y,2)*sin(2.0*M_PI*z) + 0.5*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))  ;
1143:      */
1144:     /*
1145:      Fm[0]=-2*x - 2.0*pow(M_PI,2)*sin(M_PI*x);
1146:      Fm[1]=-2*y - 2.0*pow(M_PI,2)*sin(M_PI*y);
1147:      Fm[2]=-2*z - 2.0*pow(M_PI,2)*sin(M_PI*z);
1148:      */
1149:     /*
1150:      Fm[0] = -2*x + pow(z,3)*exp(y)*sin(M_PI*x) + 6.0*z*exp(y)*sin(M_PI*x) - 6.0*M_PI*pow(x,2)*cos(2.0*M_PI*z) - 2.0*pow(M_PI,2)*pow(z,3)*exp(y)*sin(M_PI*x) + 2.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x) - 2.0*z*pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x) ;
1151:      Fm[1] = -2*y - 6.0*M_PI*pow(y,2)*cos(2.0*M_PI*z) + 2.0*z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 6.0*z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*pow(z,3)*cos(M_PI*x)*exp(y) - 4.0*M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y);
1152:      Fm[2] = -2*z - 6.0*x*sin(2.0*M_PI*z) - 6.0*y*sin(2.0*M_PI*z) - cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 8.0*pow(M_PI,2)*(pow(x,3) + pow(y,3))*sin(2.0*M_PI*z) + 3.0*M_PI*pow(z,2)*cos(M_PI*x)*exp(y) + M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;
1153:      */
1154:     Fm[0] = -2*x + 2*z*(pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x) - 1.0*M_PI*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x)) + pow(z,3)*exp(y)*sin(M_PI*x) + 6.0*z*exp(y)*sin(M_PI*x) - 1.0*pow(M_PI,2)*pow(z,3)*exp(y)*sin(M_PI*x) + 2.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x) - 2.0*z*pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x);
1155:     Fm[1] = -2*y + 2*z*(-cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)) + 2.0*z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 6.0*z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 4.0*M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;
1156:     Fm[2] = -2*z + pow(z,2)*(-2.0*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 2.0*pow(M_PI,3)*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)) + pow(z,2)*(cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 - 3*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + pow(M_PI,3)*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2 - 3*M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2) + 1.0*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 0.25*pow(M_PI,3)*pow(z,4)*cos(M_PI*x)*exp(y) - 0.25*M_PI*pow(z,4)*cos(M_PI*x)*exp(y) - 3.0*M_PI*pow(z,2)*cos(M_PI*x)*exp(y) - 1.0*M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;
1157:   }
1158:   if (Fc) {
1159:     /**Fc = -2.0*M_PI*(pow(x,3) + pow(y,3))*cos(2.0*M_PI*z) - z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*pow(z,3)*cos(M_PI*x)*exp(y) + M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;*/
1160:     /**Fc = M_PI*cos(M_PI*x) + M_PI*cos(M_PI*y) + M_PI*cos(M_PI*z);*/
1161:     /**Fc = -2.0*M_PI*(pow(x,3) + pow(y,3))*cos(2.0*M_PI*z) - z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*pow(z,3)*cos(M_PI*x)*exp(y) + M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y);*/
1162:     *Fc = 0.0;
1163:   }
1164: }

1168: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx,PetscInt my,PetscInt mz,DM *_da,Vec *_X)
1169: {
1170:   DM             da,cda;
1171:   Vec            X,local_X;
1172:   StokesDOF      ***_stokes;
1173:   Vec            coords;
1174:   DMDACoor3d     ***_coords;
1175:   PetscInt       si,sj,sk,ei,ej,ek,i,j,k;

1179:   DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1180:                       mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,4,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&da);
1181:   DMDASetFieldName(da,0,"anlytic_Vx");
1182:   DMDASetFieldName(da,1,"anlytic_Vy");
1183:   DMDASetFieldName(da,2,"anlytic_Vz");
1184:   DMDASetFieldName(da,3,"analytic_P");

1186:   DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);

1188:   DMDAGetGhostedCoordinates(da,&coords);
1189:   DMDAGetCoordinateDA(da,&cda);
1190:   DMDAVecGetArray(cda,coords,&_coords);

1192:   DMCreateGlobalVector(da,&X);
1193:   DMCreateLocalVector(da,&local_X);
1194:   DMDAVecGetArray(da,local_X,&_stokes);

1196:   DMDAGetGhostCorners(da,&si,&sj,&sk,&ei,&ej,&ek);
1197:   for (k = sk; k < sk+ek; k++) {
1198:     for (j = sj; j < sj+ej; j++) {
1199:       for (i = si; i < si+ei; i++) {
1200:         PetscReal pos[NSD],pressure,vel[NSD];

1202:         pos[0] = PetscRealPart(_coords[k][j][i].x);
1203:         pos[1] = PetscRealPart(_coords[k][j][i].y);
1204:         pos[2] = PetscRealPart(_coords[k][j][i].z);

1206:         evaluate_MS_FrankKamentski(pos,vel,&pressure,PETSC_NULL,PETSC_NULL,PETSC_NULL);

1208:         _stokes[k][j][i].u_dof = vel[0];
1209:         _stokes[k][j][i].v_dof = vel[1];
1210:         _stokes[k][j][i].w_dof = vel[2];
1211:         _stokes[k][j][i].p_dof = pressure;
1212:       }
1213:     }
1214:   }
1215:   DMDAVecRestoreArray(da,local_X,&_stokes);
1216:   DMDAVecRestoreArray(cda,coords,&_coords);

1218:   DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1219:   DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);

1221:   VecDestroy(&local_X);

1223:   *_da = da;
1224:   *_X  = X;
1225:   return(0);
1226: }

1230: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da,Vec X,Vec X_analytic)
1231: {
1232:   DM          cda;
1233:   Vec         coords,X_analytic_local,X_local;
1234:   DMDACoor3d  ***_coords;
1235:   PetscInt    sex,sey,sez,mx,my,mz;
1236:   PetscInt    ei,ej,ek;
1237:   PetscScalar el_coords[NODES_PER_EL*NSD];
1238:   StokesDOF   ***stokes_analytic,***stokes;
1239:   StokesDOF   stokes_analytic_e[NODES_PER_EL],stokes_e[NODES_PER_EL];

1241:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1242:   PetscScalar    Ni_p[NODES_PER_EL];
1243:   PetscInt       ngp;
1244:   PetscScalar    gp_xi[GAUSS_POINTS][NSD];
1245:   PetscScalar    gp_weight[GAUSS_POINTS];
1246:   PetscInt       p,i;
1247:   PetscScalar    J_p,fac;
1248:   PetscScalar    h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1249:   PetscScalar    tint_p_ms,tint_p,int_p_ms,int_p;
1250:   PetscInt       M;
1251:   PetscReal      xymin[NSD],xymax[NSD];

1255:   /* define quadrature rule */
1256:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

1258:   /* setup for coords */
1259:   DMDAGetCoordinateDA(stokes_da,&cda);
1260:   DMDAGetGhostedCoordinates(stokes_da,&coords);
1261:   DMDAVecGetArray(cda,coords,&_coords);

1263:   /* setup for analytic */
1264:   DMCreateLocalVector(stokes_da,&X_analytic_local);
1265:   DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1266:   DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1267:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1269:   /* setup for solution */
1270:   DMCreateLocalVector(stokes_da,&X_local);
1271:   DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1272:   DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1273:   DMDAVecGetArray(stokes_da,X_local,&stokes);

1275:   DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1276:   DMDAGetBoundingBox(stokes_da,xymin,xymax);

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

1280:   tp_L2 = tu_L2 = tu_H1 = 0.0;
1281:   tint_p_ms = tint_p = 0.0;

1283:   DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);

1285:   for (ek = sez; ek < sez+mz; ek++) {
1286:     for (ej = sey; ej < sey+my; ej++) {
1287:       for (ei = sex; ei < sex+mx; ei++) {
1288:         /* get coords for the element */
1289:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1290:         StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1291:         StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);

1293:         /* evaluate integral */
1294:         for (p = 0; p < ngp; p++) {
1295:           ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1296:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1297:           ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1298:           fac = gp_weight[p]*J_p;

1300:           for (i = 0; i < NODES_PER_EL; i++) {
1301:             tint_p_ms = tint_p_ms+fac*Ni_p[i]*stokes_analytic_e[i].p_dof;
1302:             tint_p    = tint_p   +fac*Ni_p[i]*stokes_e[i].p_dof;
1303:           }
1304:         }
1305:       }
1306:     }
1307:   }
1308:   MPI_Allreduce(&tint_p_ms,&int_p_ms,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1309:   MPI_Allreduce(&tint_p,&int_p,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);

1311:   PetscPrintf(PETSC_COMM_WORLD,"\\int P dv %1.4e (h)  %1.4e (ms)\n",PetscRealPart(int_p),PetscRealPart(int_p_ms));

1313:   /* remove mine and add manufacture one */
1314:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1315:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);

1317:   {
1318:     PetscInt k,L,dof;
1319:     PetscScalar *fields;
1320:     DMDAGetInfo(stokes_da,0, 0,0,0, 0,0,0, &dof,0,0,0,0,0);

1322:     VecGetLocalSize(X_local,&L);
1323:     VecGetArray(X_local,&fields);
1324:     for (k=0; k<L/dof; k++) {
1325:       fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1326:     }
1327:     VecRestoreArray(X_local,&fields);

1329:     VecGetLocalSize(X,&L);
1330:     VecGetArray(X,&fields);
1331:     for (k=0; k<L/dof; k++) {
1332:       fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1333:     }
1334:     VecRestoreArray(X,&fields);
1335:   }

1337:   DMDAVecGetArray(stokes_da,X_local,&stokes);
1338:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1340:   for (ek = sez; ek < sez+mz; ek++) {
1341:     for (ej = sey; ej < sey+my; ej++) {
1342:       for (ei = sex; ei < sex+mx; ei++) {
1343:         /* get coords for the element */
1344:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1345:         StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1346:         StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);

1348:         /* evaluate integral */
1349:         p_e_L2 = 0.0;
1350:         u_e_L2 = 0.0;
1351:         u_e_H1 = 0.0;
1352:         for (p = 0; p < ngp; p++) {
1353:           ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1354:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1355:           ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1356:           fac = gp_weight[p]*J_p;

1358:           for (i = 0; i < NODES_PER_EL; i++) {
1359:             PetscScalar u_error,v_error,w_error;

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

1363:             u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1364:             v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1365:             w_error = stokes_e[i].w_dof-stokes_analytic_e[i].w_dof;
1366:             /*
1367:              if (p==0) {
1368:              printf("p=0: %d %d %d %1.4e,%1.4e,%1.4e \n", ei,ej,ek,u_error,v_error,w_error);
1369:              }
1370:              */
1371:             u_e_L2  += fac*Ni_p[i]*(u_error*u_error+v_error*v_error+w_error*w_error);

1373:             u_e_H1 = u_e_H1+fac*( GNx_p[0][i]*u_error*GNx_p[0][i]*u_error              /* du/dx */
1374:                                   +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error
1375:                                   +GNx_p[2][i]*u_error*GNx_p[2][i]*u_error
1376:                                   +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error              /* dv/dx */
1377:                                   +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error
1378:                                   +GNx_p[2][i]*v_error*GNx_p[2][i]*v_error
1379:                                   +GNx_p[0][i]*w_error*GNx_p[0][i]*w_error              /* dw/dx */
1380:                                   +GNx_p[1][i]*w_error*GNx_p[1][i]*w_error
1381:                                   +GNx_p[2][i]*w_error*GNx_p[2][i]*w_error );
1382:           }
1383:         }

1385:         tp_L2 += p_e_L2;
1386:         tu_L2 += u_e_L2;
1387:         tu_H1 += u_e_H1;
1388:       }
1389:     }
1390:   }
1391:   MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1392:   MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1393:   MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1394:   p_L2 = PetscSqrtScalar(p_L2);
1395:   u_L2 = PetscSqrtScalar(u_L2);
1396:   u_H1 = PetscSqrtScalar(u_H1);

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


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

1403:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1404:   VecDestroy(&X_analytic_local);
1405:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1406:   VecDestroy(&X_local);
1407:   return(0);
1408: }

1412: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da,Vec FIELD,const char file_prefix[])
1413: {
1414:   char           vtk_filename[PETSC_MAX_PATH_LEN];
1415:   PetscMPIInt    rank;
1416:   MPI_Comm       comm;
1417:   FILE           *vtk_fp = PETSC_NULL;
1418:   PetscInt       si,sj,sk,nx,ny,nz,i;
1419:   PetscInt       f,n_fields,N;
1420:   DM             cda;
1421:   Vec            coords;
1422:   Vec            l_FIELD;
1423:   PetscScalar    *_L_FIELD;
1424:   PetscInt       memory_offset;
1425:   PetscScalar    *buffer;
1426:   PetscLogDouble t0,t1;

1430:   PetscGetTime(&t0);

1432:   /* create file name */
1433:   PetscObjectGetComm((PetscObject)da,&comm);
1434:   MPI_Comm_rank(comm,&rank);
1435:   PetscSNPrintf(vtk_filename,sizeof vtk_filename,"subdomain-%s-p%1.4d.vts",file_prefix,rank);

1437:   /* open file and write header */
1438:   vtk_fp = fopen(vtk_filename,"w");
1439:   if (!vtk_fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"Cannot open file = %s \n",vtk_filename);

1441:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");

1443:   /* coords */
1444:   DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1445:   N = nx * ny * nz;

1447: #ifdef PETSC_WORDS_BIGENDIAN
1448:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1449: #else
1450:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1451: #endif
1452:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",si,si+nx-1,sj,sj+ny-1,sk,sk+nz-1);
1453:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    <Piece Extent=\"%D %D %D %D %D %D\">\n",si,si+nx-1,sj,sj+ny-1,sk,sk+nz-1);

1455:   memory_offset = 0;

1457:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      <CellData></CellData>\n");

1459:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      <Points>\n");

1461:   /* copy coordinates */
1462:   DMDAGetCoordinateDA(da,&cda);
1463:   DMDAGetGhostedCoordinates(da,&coords);
1464:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"        <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%d\" />\n",memory_offset);
1465:   memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N*3;

1467:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      </Points>\n");

1469:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      <PointData Scalars=\" ");
1470:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_fields,0,0,0,0,0);
1471:   for (f=0; f<n_fields; f++) {
1472:     const char *field_name;
1473:     DMDAGetFieldName(da,f,&field_name);
1474:     PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"%s ",field_name);
1475:   }
1476:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\">\n");

1478:   for (f=0; f<n_fields; f++) {
1479:     const char *field_name;

1481:     DMDAGetFieldName(da,f,&field_name);
1482:     PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"        <DataArray type=\"Float64\" Name=\"%s\" format=\"appended\" offset=\"%d\"/>\n", field_name,memory_offset);
1483:     memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N;
1484:   }

1486:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      </PointData>\n");
1487:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </Piece>\n");
1488:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  </StructuredGrid>\n");

1490:   PetscMalloc(sizeof(PetscScalar)*N,&buffer);
1491:   DMGetLocalVector(da,&l_FIELD);
1492:   DMGlobalToLocalBegin(da, FIELD,INSERT_VALUES,l_FIELD);
1493:   DMGlobalToLocalEnd(da,FIELD,INSERT_VALUES,l_FIELD);
1494:   VecGetArray(l_FIELD,&_L_FIELD);

1496:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  <AppendedData encoding=\"raw\">\n");
1497:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"_");

1499:   /* write coordinates */
1500:   {
1501:     int length = sizeof(PetscScalar)*N*3;
1502:     PetscScalar *allcoords;

1504:     fwrite(&length,sizeof(int),1,vtk_fp);
1505:     VecGetArray(coords,&allcoords);
1506:     fwrite(allcoords,sizeof(PetscScalar),3*N,vtk_fp);
1507:     VecRestoreArray(coords,&allcoords);
1508:   }
1509:   /* write fields */
1510:   for (f=0; f<n_fields; f++) {
1511:     int length = sizeof(PetscScalar)*N;
1512:     fwrite(&length,sizeof(int),1,vtk_fp);
1513:     /* load */
1514:     for (i=0; i<N; i++) {
1515:       buffer[i] = _L_FIELD[n_fields*i + f];
1516:     }

1518:     /* write */
1519:     fwrite(buffer,sizeof(PetscScalar),N,vtk_fp);
1520:   }
1521:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\n  </AppendedData>\n");

1523:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");

1525:   PetscFree(buffer);
1526:   VecRestoreArray(l_FIELD,&_L_FIELD);
1527:   DMRestoreLocalVector(da,&l_FIELD);

1529:   if (vtk_fp) {
1530:     fclose(vtk_fp);
1531:     vtk_fp = NULL;
1532:   }

1534:   PetscGetTime(&t1);
1535:   return(0);
1536: }

1540: PetscErrorCode DAViewVTK_write_PieceExtend(FILE *vtk_fp,PetscInt indent_level,DM da,const char local_file_prefix[])
1541: {
1542:   PetscMPIInt    nproc,rank;
1543:   MPI_Comm       comm;
1544:   const PetscInt *lx,*ly,*lz;
1545:   PetscInt       M,N,P,pM,pN,pP,sum,*olx,*oly,*olz;
1546:   PetscInt       *osx,*osy,*osz,*oex,*oey,*oez;
1547:   PetscInt       i,j,k,II,stencil;

1551:   /* create file name */
1552:   PetscObjectGetComm((PetscObject)da,&comm);
1553:   MPI_Comm_size(comm,&nproc);
1554:   MPI_Comm_rank(comm,&rank);

1556:   DMDAGetInfo(da,0,&M,&N,&P,&pM,&pN,&pP,0,&stencil,0,0,0,0);
1557:   DMDAGetOwnershipRanges(da,&lx,&ly,&lz);

1559:   /* generate start,end list */
1560:   PetscMalloc(sizeof(PetscInt)*(pM+1),&olx);
1561:   PetscMalloc(sizeof(PetscInt)*(pN+1),&oly);
1562:   PetscMalloc(sizeof(PetscInt)*(pP+1),&olz);
1563:   sum = 0;
1564:   for (i=0; i<pM; i++) {
1565:     olx[i] = sum;
1566:     sum = sum + lx[i];
1567:   } olx[pM] = sum;
1568:   sum = 0;
1569:   for (i=0; i<pN; i++) {
1570:     oly[i] = sum;
1571:     sum = sum + ly[i];
1572:   } oly[pN] = sum;
1573:   sum = 0;
1574:   for (i=0; i<pP; i++) {
1575:     olz[i] = sum;
1576:     sum = sum + lz[i];
1577:   } olz[pP] = sum;
1578:   PetscMalloc(sizeof(PetscInt)*(pM),&osx);
1579:   PetscMalloc(sizeof(PetscInt)*(pN),&osy);
1580:   PetscMalloc(sizeof(PetscInt)*(pP),&osz);
1581:   PetscMalloc(sizeof(PetscInt)*(pM),&oex);
1582:   PetscMalloc(sizeof(PetscInt)*(pN),&oey);
1583:   PetscMalloc(sizeof(PetscInt)*(pP),&oez);
1584:   for (i=0; i<pM; i++) {
1585:     osx[i] = olx[i] - stencil;
1586:     oex[i] = olx[i] + lx[i] + stencil;
1587:     if (osx[i]<0)osx[i]=0;
1588:     if (oex[i]>M)oex[i]=M;
1589:   }

1591:   for (i=0; i<pN; i++) {
1592:     osy[i] = oly[i] - stencil;
1593:     oey[i] = oly[i] + ly[i] + stencil;
1594:     if (osy[i]<0)osy[i]=0;
1595:     if (oey[i]>M)oey[i]=N;
1596:   }
1597:   for (i=0; i<pP; i++) {
1598:     osz[i] = olz[i] - stencil;
1599:     oez[i] = olz[i] + lz[i] + stencil;
1600:     if (osz[i]<0)osz[i]=0;
1601:     if (oez[i]>P)oez[i]=P;
1602:   }

1604:   for (k=0; k<pP; k++) {
1605:     for (j=0; j<pN; j++) {
1606:       for (i=0; i<pM; i++) {
1607:         char name[PETSC_MAX_PATH_LEN];
1608:         PetscInt procid = i + j*pM + k*pM*pN; /* convert proc(i,j,k) to pid */
1609:         PetscSNPrintf(name,sizeof name,"subdomain-%s-p%1.4d.vts",local_file_prefix,procid);
1610:         for (II=0; II<indent_level; II++) {
1611:           PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  ");
1612:         }
1613:         PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<Piece Extent=\"%d %d %d %d %d %d\"      Source=\"%s\"/>\n",
1614:                      osx[i],oex[i]-1,
1615:                      osy[j],oey[j]-1,
1616:                      osz[k],oez[k]-1,name);
1617:       }
1618:     }
1619:   }
1620:   PetscFree(olx);
1621:   PetscFree(oly);
1622:   PetscFree(olz);
1623:   PetscFree(osx);
1624:   PetscFree(osy);
1625:   PetscFree(osz);
1626:   PetscFree(oex);
1627:   PetscFree(oey);
1628:   PetscFree(oez);
1629:   return(0);
1630: }

1634: PetscErrorCode DAView_3DVTK_PStructuredGrid(DM da,const char file_prefix[],const char local_file_prefix[])
1635: {
1636:   MPI_Comm       comm;
1637:   PetscMPIInt    nproc,rank;
1638:   char           vtk_filename[PETSC_MAX_PATH_LEN];
1639:   FILE           *vtk_fp = PETSC_NULL;
1640:   PetscInt       M,N,P,si,sj,sk,nx,ny,nz;
1641:   PetscInt       i,dofs;

1645:   /* only master generates this file */
1646:   PetscObjectGetComm((PetscObject)da,&comm);
1647:   MPI_Comm_size(comm,&nproc);
1648:   MPI_Comm_rank(comm,&rank);

1650:   if (rank != 0) { return(0); }

1652:   /* create file name */
1653:   PetscSNPrintf(vtk_filename,sizeof vtk_filename,"%s.pvts",file_prefix);
1654:   vtk_fp = fopen(vtk_filename,"w");
1655:   if (!vtk_fp) { SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"Cannot open file = %s \n",vtk_filename); }

1657:   /* (VTK) generate pvts header */
1658:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");

1660: #ifdef PETSC_WORDS_BIGENDIAN
1661:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1662: #else
1663:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1664: #endif

1666:   /* define size of the nodal mesh based on the cell DM */
1667:   DMDAGetInfo(da,0,&M,&N,&P,0,0,0,&dofs,0,0,0,0,0);
1668:   DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1669:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  <PStructuredGrid GhostLevel=\"1\" WholeExtent=\"%d %d %d %d %d %d\">\n",0,M-1,0,N-1,0,P-1); /* note overlap = 1 for Q1 */

1671:   /* DUMP THE CELL REFERENCES */
1672:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    <PCellData>\n");
1673:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </PCellData>\n");

1675:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    <PPoints>\n");
1676:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      <PDataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"%d\"/>\n",NSD);
1677:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </PPoints>\n");

1679:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    <PPointData>\n");
1680:   for (i=0; i<dofs; i++) {
1681:     const char *fieldname;
1682:     DMDAGetFieldName(da,i,&fieldname);
1683:     PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n",fieldname);
1684:   }
1685:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </PPointData>\n");

1687:   /* write out the parallel information */
1688:   DAViewVTK_write_PieceExtend(vtk_fp,2,da,local_file_prefix);

1690:   /* close the file */
1691:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  </PStructuredGrid>\n");
1692:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");

1694:   if (vtk_fp){
1695:     fclose(vtk_fp);
1696:     vtk_fp = NULL;
1697:   }
1698:   return(0);
1699: }

1703: PetscErrorCode DAView3DPVTS(DM da, Vec x,const char NAME[])
1704: {
1705:   char           vts_filename[PETSC_MAX_PATH_LEN];
1706:   char           pvts_filename[PETSC_MAX_PATH_LEN];

1710:   PetscSNPrintf(vts_filename,sizeof vts_filename,"%s-mesh",NAME);
1711:   DAView_3DVTK_StructuredGrid_appended(da,x,vts_filename);

1713:   PetscSNPrintf(pvts_filename,sizeof pvts_filename,"%s-mesh",NAME);
1714:   DAView_3DVTK_PStructuredGrid(da,pvts_filename,vts_filename);
1715:   return(0);
1716: }

1720: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
1721: {
1723:   PetscReal      norms[4];
1724:   Vec            Br,v,w;
1725:   Mat            A;

1728:   KSPGetOperators(ksp,&A,0,0);
1729:   MatGetVecs(A,&w,&v);

1731:   KSPBuildResidual(ksp,v,w,&Br);

1733:   VecStrideNorm(Br,0,NORM_2,&norms[0]);
1734:   VecStrideNorm(Br,1,NORM_2,&norms[1]);
1735:   VecStrideNorm(Br,2,NORM_2,&norms[2]);
1736:   VecStrideNorm(Br,3,NORM_2,&norms[3]);

1738:   VecDestroy(&v);
1739:   VecDestroy(&w);

1741:   PetscPrintf(PETSC_COMM_WORLD,"%3D KSP Component U,V,W,P residual norm [ %1.12e, %1.12e, %1.12e, %1.12e ]\n",n,norms[0],norms[1],norms[2],norms[3]);
1742:   return(0);
1743: }

1747: static PetscErrorCode PCMGSetupViaCoarsen(PC pc,DM da_fine)
1748: {
1749:   PetscInt       nlevels,k,PETSC_UNUSED finest;
1750:   DM             *da_list,*daclist;
1751:   Mat            R;

1755:   nlevels = 1;
1756:   PetscOptionsGetInt(PETSC_NULL,"-levels",&nlevels,0);

1758:   PetscMalloc(sizeof(DM)*nlevels,&da_list);
1759:   for (k=0; k<nlevels; k++) {  da_list[k] = PETSC_NULL;  }
1760:   PetscMalloc(sizeof(DM)*nlevels,&daclist);
1761:   for (k=0; k<nlevels; k++) {  daclist[k] = PETSC_NULL;  }

1763:   /* finest grid is nlevels - 1 */
1764:   finest = nlevels - 1;
1765:   daclist[0] = da_fine;
1766:   PetscObjectReference((PetscObject)da_fine);
1767:   DMCoarsenHierarchy(da_fine,nlevels-1,&daclist[1]);
1768:   for (k=0; k<nlevels; k++) {
1769:     da_list[k] = daclist[nlevels-1-k];
1770:     DMDASetUniformCoordinates(da_list[k],0.0,1.0,0.0,1.0,0.0,1.0);
1771:   }

1773:   PCMGSetLevels(pc,nlevels,PETSC_NULL);
1774:   PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
1775:   PCMGSetGalerkin(pc,PETSC_TRUE);

1777:   for (k=1; k<nlevels; k++){
1778:     DMCreateInterpolation(da_list[k-1],da_list[k],&R,PETSC_NULL);
1779:     PCMGSetInterpolation(pc,k,R);
1780:     MatDestroy(&R);
1781:   }

1783:   /* tidy up */
1784:   for (k=0; k<nlevels; k++) {
1785:     DMDestroy(&da_list[k]);
1786:   }
1787:   PetscFree(da_list);
1788:   PetscFree(daclist);
1789:   return(0);
1790: }

1794: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx,PetscInt my,PetscInt mz)
1795: {
1796:   DM                     da_Stokes;
1797:   PetscInt               u_dof,p_dof,dof,stencil_width;
1798:   Mat                    A,B;
1799:   PetscInt               mxl,myl,mzl;
1800:   DM                     vel_cda;
1801:   Vec                    vel_coords;
1802:   PetscInt               p;
1803:   Vec                    f,X;
1804:   DMDACoor3d             ***_vel_coords;
1805:   PetscInt               its;
1806:   KSP                    ksp_S;
1807:   PetscInt               model_definition = 0;
1808:   PetscInt               ei,ej,ek,sex,sey,sez,Imx,Jmy,Kmz;
1809:   CellProperties         cell_properties;
1810:   PetscBool              write_output = PETSC_FALSE;
1811:   PetscErrorCode         ierr;

1814:   /* Generate the da for velocity and pressure */
1815:   /* Num nodes in each direction is mx+1, my+1, mz+1 */
1816:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1817:   p_dof         = P_DOFS; /* p - pressure */
1818:   dof           = u_dof+p_dof;
1819:   stencil_width = 1;
1820:   DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1821:                                mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,PETSC_NULL,PETSC_NULL,PETSC_NULL,&da_Stokes);
1822:   DMDASetFieldName(da_Stokes,0,"Vx");
1823:   DMDASetFieldName(da_Stokes,1,"Vy");
1824:   DMDASetFieldName(da_Stokes,2,"Vz");
1825:   DMDASetFieldName(da_Stokes,3,"P");

1827:   /* unit box [0,1] x [0,1] x [0,1] */
1828:   DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.0,1.0);

1830:   /* local number of elements */
1831:   DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,&mzl);

1833:   /* create quadrature point info for PDE coefficients */
1834:   CellPropertiesCreate(da_Stokes,&cell_properties);

1836:   /* interpolate the coordinates to quadrature points */
1837:   DMDAGetCoordinateDA(da_Stokes,&vel_cda);
1838:   DMDAGetGhostedCoordinates(da_Stokes,&vel_coords);
1839:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1840:   DMDAGetElementCorners(da_Stokes,&sex,&sey,&sez,&Imx,&Jmy,&Kmz);
1841:   for (ek = sez; ek < sez+Kmz; ek++) {
1842:     for (ej = sey; ej < sey+Jmy; ej++) {
1843:       for (ei = sex; ei < sex+Imx; ei++) {
1844:         /* get coords for the element */
1845:         PetscInt    ngp;
1846:         PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
1847:         PetscScalar el_coords[NSD*NODES_PER_EL];
1848:         GaussPointCoefficients *cell;

1850:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1851:         GetElementCoords3D(_vel_coords,ei,ej,ek,el_coords);
1852:         ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

1854:         for (p = 0; p < GAUSS_POINTS; p++) {
1855:           PetscScalar xi_p[NSD],Ni_p[NODES_PER_EL];
1856:           PetscScalar gp_x,gp_y,gp_z;
1857:           PetscInt    n;

1859:           xi_p[0] = gp_xi[p][0];
1860:           xi_p[1] = gp_xi[p][1];
1861:           xi_p[2] = gp_xi[p][2];
1862:           ShapeFunctionQ13D_Evaluate(xi_p,Ni_p);

1864:           gp_x = gp_y = gp_z = 0.0;
1865:           for (n = 0; n < NODES_PER_EL; n++) {
1866:             gp_x = gp_x+Ni_p[n]*el_coords[NSD*n  ];
1867:             gp_y = gp_y+Ni_p[n]*el_coords[NSD*n+1];
1868:             gp_z = gp_z+Ni_p[n]*el_coords[NSD*n+2];
1869:           }
1870:           cell->gp_coords[NSD*p  ] = gp_x;
1871:           cell->gp_coords[NSD*p+1] = gp_y;
1872:           cell->gp_coords[NSD*p+2] = gp_z;
1873:         }
1874:       }
1875:     }
1876:   }

1878:   PetscOptionsGetInt(PETSC_NULL,"-model",&model_definition,PETSC_NULL);

1880:   switch (model_definition) {
1881:   case 0: /* isoviscous */
1882:     for (ek = sez; ek < sez+Kmz; ek++) {
1883:       for (ej = sey; ej < sey+Jmy; ej++) {
1884:         for (ei = sex; ei < sex+Imx; ei++) {
1885:           GaussPointCoefficients *cell;

1887:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1888:           for (p = 0; p < GAUSS_POINTS; p++) {
1889:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p  ]);
1890:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1891:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);

1893:             cell->eta[p] = 1.0;

1895:             cell->fx[p] = 0.0*coord_x;
1896:             cell->fy[p] = -sin((double)2.2*M_PI*coord_y)*cos(1.0*M_PI*coord_x);
1897:             cell->fz[p] = 0.0*coord_z;
1898:             cell->hc[p] = 0.0;
1899:           }
1900:         }
1901:       }
1902:     }
1903:     break;

1905:   case 1: /* manufactured */
1906:     for (ek = sez; ek < sez+Kmz; ek++) {
1907:       for (ej = sey; ej < sey+Jmy; ej++) {
1908:         for (ei = sex; ei < sex+Imx; ei++) {
1909:           PetscReal eta,Fm[NSD],Fc,pos[NSD];
1910:           GaussPointCoefficients *cell;

1912:           CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1913:           for (p = 0; p < GAUSS_POINTS; p++) {
1914:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p  ]);
1915:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1916:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);

1918:             pos[0] = coord_x;
1919:             pos[1] = coord_y;
1920:             pos[2] = coord_z;

1922:             evaluate_MS_FrankKamentski(pos,PETSC_NULL,PETSC_NULL,&eta,Fm,&Fc);
1923:             cell->eta[p] = eta;

1925:             cell->fx[p] = Fm[0];
1926:             cell->fy[p] = Fm[1];
1927:             cell->fz[p] = Fm[2];
1928:             cell->hc[p] = Fc;
1929:           }
1930:         }
1931:       }
1932:     }
1933:     break;

1935:   case 2: /* solcx */
1936:     for (ek = sez; ek < sez+Kmz; ek++) {
1937:       for (ej = sey; ej < sey+Jmy; ej++) {
1938:         for (ei = sex; ei < sex+Imx; ei++) {
1939:           GaussPointCoefficients *cell;

1941:           CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1942:           for (p = 0; p < GAUSS_POINTS; p++) {
1943:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p  ]);
1944:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1945:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);

1947:             cell->eta[p] = 1.0;

1949:             cell->fx[p] = 0.0;
1950:             cell->fy[p] = -sin((double)3*M_PI*coord_y)*cos(1.0*M_PI*coord_x);
1951:             cell->fz[p] = 0.0*coord_z;
1952:             cell->hc[p] = 0.0;
1953:           }
1954:         }
1955:       }
1956:     }
1957:     break;

1959:   case 3: /* sinker */
1960:     for (ek = sez; ek < sez+Kmz; ek++) {
1961:       for (ej = sey; ej < sey+Jmy; ej++) {
1962:         for (ei = sex; ei < sex+Imx; ei++) {
1963:           GaussPointCoefficients *cell;

1965:           CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1966:           for (p = 0; p < GAUSS_POINTS; p++) {
1967:             PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p  ]);
1968:             PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);
1969:             PetscReal zp = PetscRealPart(cell->gp_coords[NSD*p+2]);

1971:             cell->eta[p] = 1.0e-2;
1972:             cell->fx[p] = 0.0;
1973:             cell->fy[p] = 0.0;
1974:             cell->fz[p] = 0.0;
1975:             cell->hc[p] = 0.0;

1977:             if ( (fabs(xp-0.5) < 0.2) &&
1978:                  (fabs(yp-0.5) < 0.2) &&
1979:                  (fabs(zp-0.5) < 0.2) ) {
1980:               cell->eta[p] = 1.0;
1981:               cell->fz[p]  = 1.0;
1982:             }

1984:           }
1985:         }
1986:       }
1987:     }
1988:     break;

1990:   default:
1991:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"No default model is supported. Choose either -model {0,1,2,3}");
1992:     break;
1993:   }

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

1997:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1998:   DMCreateMatrix(da_Stokes,MATAIJ,&A);
1999:   DMCreateMatrix(da_Stokes,MATAIJ,&B);
2000:   DMCreateGlobalVector(da_Stokes,&X);
2001:   DMCreateGlobalVector(da_Stokes,&f);

2003:   /* assemble A11 */
2004:   MatZeroEntries(A);
2005:   MatZeroEntries(B);
2006:   VecZeroEntries(f);

2008:   AssembleA_Stokes(A,da_Stokes,cell_properties);
2009:   AssembleA_PCStokes(B,da_Stokes,cell_properties);
2010:   /* build force vector */
2011:   AssembleF_Stokes(f,da_Stokes,cell_properties);

2013:   /* SOLVE */
2014:   KSPCreate(PETSC_COMM_WORLD,&ksp_S);
2015:   KSPSetOptionsPrefix(ksp_S,"stokes_"); /* stokes */
2016:   KSPSetOperators(ksp_S,A,B,SAME_NONZERO_PATTERN);
2017:   KSPSetFromOptions(ksp_S);

2019:   {
2020:     PC pc;
2021:     const PetscInt ufields[] = {0,1,2},pfields[] = {3};
2022:     KSPGetPC(ksp_S,&pc);
2023:     PCFieldSplitSetBlockSize(pc,4);
2024:     PCFieldSplitSetFields(pc,"u",3,ufields,ufields);
2025:     PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
2026:   }

2028:   {
2029:     PC pc;
2030:     PetscBool same = PETSC_FALSE;
2031:     KSPGetPC(ksp_S,&pc);
2032:     PetscObjectTypeCompare((PetscObject)pc,PCMG,&same);
2033:     if (same) {
2034:       PCMGSetupViaCoarsen(pc,da_Stokes);
2035:     }
2036:   }

2038:   {
2039:     PetscBool stokes_monitor = PETSC_FALSE;
2040:     PetscOptionsGetBool(PETSC_NULL,"-stokes_ksp_monitor_blocks",&stokes_monitor,0);
2041:     if (stokes_monitor) {
2042:       KSPMonitorSet(ksp_S,KSPMonitorStokesBlocks,PETSC_NULL,PETSC_NULL);
2043:     }
2044:   }
2045:   KSPSolve(ksp_S,f,X);

2047:   PetscOptionsGetBool(PETSC_NULL,"-write_pvts",&write_output,PETSC_NULL);
2048:   if (write_output) {DAView3DPVTS(da_Stokes,X,"up");}
2049:   {
2050:     PetscBool flg = PETSC_FALSE;
2051:     char filename[PETSC_MAX_PATH_LEN];
2052:     PetscOptionsGetString(PETSC_NULL,"-write_binary",filename,sizeof filename,&flg);
2053:     if (flg) {
2054:       PetscViewer viewer;
2055:       //PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer);
2056:       PetscViewerVTKOpen(PETSC_COMM_WORLD,"ex42.vts",FILE_MODE_WRITE,&viewer);
2057:       VecView(X,viewer);
2058:       PetscViewerDestroy(&viewer);
2059:     }
2060:   }
2061:   KSPGetIterationNumber(ksp_S,&its);

2063:   /* verify */
2064:   if (model_definition == 1) {
2065:     DM  da_Stokes_analytic;
2066:     Vec X_analytic;

2068:     DMDACreateManufacturedSolution(mx,my,mz,&da_Stokes_analytic,&X_analytic);
2069:     if (write_output) {
2070:       DAView3DPVTS(da_Stokes_analytic,X_analytic,"ms");
2071:     }
2072:     DMDAIntegrateErrors3D(da_Stokes_analytic,X,X_analytic);
2073:     if (write_output) {
2074:       DAView3DPVTS(da_Stokes,X,"up2");
2075:     }
2076:     DMDestroy(&da_Stokes_analytic);
2077:     VecDestroy(&X_analytic);
2078:   }

2080:   KSPDestroy(&ksp_S);
2081:   VecDestroy(&X);
2082:   VecDestroy(&f);
2083:   MatDestroy(&A);
2084:   MatDestroy(&B);

2086:   CellPropertiesDestroy(&cell_properties);
2087:   DMDestroy(&da_Stokes);
2088:   return(0);
2089: }

2093: int main(int argc,char **args)
2094: {
2096:   PetscInt       mx,my,mz;

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

2100:   mx = my = mz = 10;
2101:   PetscOptionsGetInt(PETSC_NULL,"-mx",&mx,PETSC_NULL);
2102:   my = mx; mz = mx;
2103:   PetscOptionsGetInt(PETSC_NULL,"-my",&my,PETSC_NULL);
2104:   PetscOptionsGetInt(PETSC_NULL,"-mz",&mz,PETSC_NULL);

2106:   solve_stokes_3d_coupled(mx,my,mz);

2108:   PetscFinalize();
2109:   return 0;
2110: }