Actual source code: ex42.c

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

 10: /* Contributed by Dave May */

 12: #include <petscksp.h>
 13: #include <petscdmda.h>

 15: #define PROFILE_TIMING
 16: #define ASSEMBLE_LOWER_TRIANGULAR

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

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

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

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

 49: /* elements */
 50: PetscErrorCode CellPropertiesCreate(DM da_stokes, CellProperties *C)
 51: {
 52:   CellProperties cells;
 53:   PetscInt       mx, my, mz, sex, sey, sez;

 56:   PetscNew(&cells);

 58:   DMDAGetElementsCorners(da_stokes, &sex, &sey, &sez);
 59:   DMDAGetElementsSizes(da_stokes, &mx, &my, &mz);
 60:   cells->mx     = mx;
 61:   cells->my     = my;
 62:   cells->mz     = mz;
 63:   cells->ncells = mx * my * mz;
 64:   cells->sex    = sex;
 65:   cells->sey    = sey;
 66:   cells->sez    = sez;

 68:   PetscMalloc1(mx * my * mz, &cells->gpc);

 70:   *C = cells;
 71:   return 0;
 72: }

 74: PetscErrorCode CellPropertiesDestroy(CellProperties *C)
 75: {
 76:   CellProperties cells;

 79:   if (!C) return 0;
 80:   cells = *C;
 81:   PetscFree(cells->gpc);
 82:   PetscFree(cells);
 83:   *C = NULL;
 84:   return 0;
 85: }

 87: PetscErrorCode CellPropertiesGetCell(CellProperties C, PetscInt II, PetscInt J, PetscInt K, GaussPointCoefficients **G)
 88: {
 90:   *G = &C->gpc[(II - C->sex) + (J - C->sey) * C->mx + (K - C->sez) * C->mx * C->my];
 91:   return 0;
 92: }

 94: /* FEM routines */
 95: /*
 96:  Element: Local basis function ordering
 97:  1-----2
 98:  |     |
 99:  |     |
100:  0-----3
101:  */
102: static void ShapeFunctionQ13D_Evaluate(PetscScalar _xi[], PetscScalar Ni[])
103: {
104:   PetscReal xi   = PetscRealPart(_xi[0]);
105:   PetscReal eta  = PetscRealPart(_xi[1]);
106:   PetscReal zeta = PetscRealPart(_xi[2]);

108:   Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
109:   Ni[1] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
110:   Ni[2] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
111:   Ni[3] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);

113:   Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
114:   Ni[5] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
115:   Ni[6] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
116:   Ni[7] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
117: }

119: static void ShapeFunctionQ13D_Evaluate_dxi(PetscScalar _xi[], PetscScalar GNi[][NODES_PER_EL])
120: {
121:   PetscReal xi   = PetscRealPart(_xi[0]);
122:   PetscReal eta  = PetscRealPart(_xi[1]);
123:   PetscReal zeta = PetscRealPart(_xi[2]);
124:   /* xi deriv */
125:   GNi[0][0] = -0.125 * (1.0 - eta) * (1.0 - zeta);
126:   GNi[0][1] = -0.125 * (1.0 + eta) * (1.0 - zeta);
127:   GNi[0][2] = 0.125 * (1.0 + eta) * (1.0 - zeta);
128:   GNi[0][3] = 0.125 * (1.0 - eta) * (1.0 - zeta);

130:   GNi[0][4] = -0.125 * (1.0 - eta) * (1.0 + zeta);
131:   GNi[0][5] = -0.125 * (1.0 + eta) * (1.0 + zeta);
132:   GNi[0][6] = 0.125 * (1.0 + eta) * (1.0 + zeta);
133:   GNi[0][7] = 0.125 * (1.0 - eta) * (1.0 + zeta);
134:   /* eta deriv */
135:   GNi[1][0] = -0.125 * (1.0 - xi) * (1.0 - zeta);
136:   GNi[1][1] = 0.125 * (1.0 - xi) * (1.0 - zeta);
137:   GNi[1][2] = 0.125 * (1.0 + xi) * (1.0 - zeta);
138:   GNi[1][3] = -0.125 * (1.0 + xi) * (1.0 - zeta);

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

150:   GNi[2][4] = 0.125 * (1.0 - xi) * (1.0 - eta);
151:   GNi[2][5] = 0.125 * (1.0 - xi) * (1.0 + eta);
152:   GNi[2][6] = 0.125 * (1.0 + xi) * (1.0 + eta);
153:   GNi[2][7] = 0.125 * (1.0 + xi) * (1.0 - eta);
154: }

156: static void matrix_inverse_3x3(PetscScalar A[3][3], PetscScalar B[3][3])
157: {
158:   PetscScalar t4, t6, t8, t10, t12, t14, t17;

160:   t4  = A[2][0] * A[0][1];
161:   t6  = A[2][0] * A[0][2];
162:   t8  = A[1][0] * A[0][1];
163:   t10 = A[1][0] * A[0][2];
164:   t12 = A[0][0] * A[1][1];
165:   t14 = A[0][0] * A[1][2];
166:   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]);

168:   B[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) * t17;
169:   B[0][1] = -(A[0][1] * A[2][2] - A[0][2] * A[2][1]) * t17;
170:   B[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) * t17;
171:   B[1][0] = -(-A[2][0] * A[1][2] + A[1][0] * A[2][2]) * t17;
172:   B[1][1] = (-t6 + A[0][0] * A[2][2]) * t17;
173:   B[1][2] = -(-t10 + t14) * t17;
174:   B[2][0] = (-A[2][0] * A[1][1] + A[1][0] * A[2][1]) * t17;
175:   B[2][1] = -(-t4 + A[0][0] * A[2][1]) * t17;
176:   B[2][2] = (-t8 + t12) * t17;
177: }

179: static void ShapeFunctionQ13D_Evaluate_dx(PetscScalar GNi[][NODES_PER_EL], PetscScalar GNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
180: {
181:   PetscScalar J00, J01, J02, J10, J11, J12, J20, J21, J22;
182:   PetscInt    n;
183:   PetscScalar iJ[3][3], JJ[3][3];

185:   J00 = J01 = J02 = 0.0;
186:   J10 = J11 = J12 = 0.0;
187:   J20 = J21 = J22 = 0.0;
188:   for (n = 0; n < NODES_PER_EL; n++) {
189:     PetscScalar cx = coords[NSD * n + 0];
190:     PetscScalar cy = coords[NSD * n + 1];
191:     PetscScalar cz = coords[NSD * n + 2];

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

198:     J10 = J10 + GNi[1][n] * cx; /* J_yx = dy/dxi */
199:     J11 = J11 + GNi[1][n] * cy; /* J_yy */
200:     J12 = J12 + GNi[1][n] * cz; /* J_yz */

202:     J20 = J20 + GNi[2][n] * cx; /* J_zx */
203:     J21 = J21 + GNi[2][n] * cy; /* J_zy */
204:     J22 = J22 + GNi[2][n] * cz; /* J_zz */
205:   }

207:   JJ[0][0] = J00;
208:   JJ[0][1] = J01;
209:   JJ[0][2] = J02;
210:   JJ[1][0] = J10;
211:   JJ[1][1] = J11;
212:   JJ[1][2] = J12;
213:   JJ[2][0] = J20;
214:   JJ[2][1] = J21;
215:   JJ[2][2] = J22;

217:   matrix_inverse_3x3(JJ, iJ);

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

221:   for (n = 0; n < NODES_PER_EL; n++) {
222:     GNx[0][n] = GNi[0][n] * iJ[0][0] + GNi[1][n] * iJ[0][1] + GNi[2][n] * iJ[0][2];
223:     GNx[1][n] = GNi[0][n] * iJ[1][0] + GNi[1][n] * iJ[1][1] + GNi[2][n] * iJ[1][2];
224:     GNx[2][n] = GNi[0][n] * iJ[2][0] + GNi[1][n] * iJ[2][1] + GNi[2][n] * iJ[2][2];
225:   }
226: }

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

244:   gp_xi[4][0] = -0.57735026919;
245:   gp_xi[4][1] = -0.57735026919;
246:   gp_xi[4][2] = 0.57735026919;
247:   gp_xi[5][0] = -0.57735026919;
248:   gp_xi[5][1] = 0.57735026919;
249:   gp_xi[5][2] = 0.57735026919;
250:   gp_xi[6][0] = 0.57735026919;
251:   gp_xi[6][1] = 0.57735026919;
252:   gp_xi[6][2] = 0.57735026919;
253:   gp_xi[7][0] = 0.57735026919;
254:   gp_xi[7][1] = -0.57735026919;
255:   gp_xi[7][2] = 0.57735026919;

257:   gp_weight[0] = 1.0;
258:   gp_weight[1] = 1.0;
259:   gp_weight[2] = 1.0;
260:   gp_weight[3] = 1.0;

262:   gp_weight[4] = 1.0;
263:   gp_weight[5] = 1.0;
264:   gp_weight[6] = 1.0;
265:   gp_weight[7] = 1.0;
266: }

268: /*
269:  i,j are the element indices
270:  The unknown is a vector quantity.
271:  The s[].c is used to indicate the degree of freedom.
272:  */
273: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[], MatStencil s_p[], PetscInt i, PetscInt j, PetscInt k)
274: {
275:   PetscInt n;

278:   /* velocity */
279:   n = 0;
280:   /* node 0 */
281:   s_u[n].i = i;
282:   s_u[n].j = j;
283:   s_u[n].k = k;
284:   s_u[n].c = 0;
285:   n++; /* Vx0 */
286:   s_u[n].i = i;
287:   s_u[n].j = j;
288:   s_u[n].k = k;
289:   s_u[n].c = 1;
290:   n++; /* Vy0 */
291:   s_u[n].i = i;
292:   s_u[n].j = j;
293:   s_u[n].k = k;
294:   s_u[n].c = 2;
295:   n++; /* Vz0 */

297:   s_u[n].i = i;
298:   s_u[n].j = j + 1;
299:   s_u[n].k = k;
300:   s_u[n].c = 0;
301:   n++;
302:   s_u[n].i = i;
303:   s_u[n].j = j + 1;
304:   s_u[n].k = k;
305:   s_u[n].c = 1;
306:   n++;
307:   s_u[n].i = i;
308:   s_u[n].j = j + 1;
309:   s_u[n].k = k;
310:   s_u[n].c = 2;
311:   n++;

313:   s_u[n].i = i + 1;
314:   s_u[n].j = j + 1;
315:   s_u[n].k = k;
316:   s_u[n].c = 0;
317:   n++;
318:   s_u[n].i = i + 1;
319:   s_u[n].j = j + 1;
320:   s_u[n].k = k;
321:   s_u[n].c = 1;
322:   n++;
323:   s_u[n].i = i + 1;
324:   s_u[n].j = j + 1;
325:   s_u[n].k = k;
326:   s_u[n].c = 2;
327:   n++;

329:   s_u[n].i = i + 1;
330:   s_u[n].j = j;
331:   s_u[n].k = k;
332:   s_u[n].c = 0;
333:   n++;
334:   s_u[n].i = i + 1;
335:   s_u[n].j = j;
336:   s_u[n].k = k;
337:   s_u[n].c = 1;
338:   n++;
339:   s_u[n].i = i + 1;
340:   s_u[n].j = j;
341:   s_u[n].k = k;
342:   s_u[n].c = 2;
343:   n++;

345:   /* */
346:   s_u[n].i = i;
347:   s_u[n].j = j;
348:   s_u[n].k = k + 1;
349:   s_u[n].c = 0;
350:   n++; /* Vx4 */
351:   s_u[n].i = i;
352:   s_u[n].j = j;
353:   s_u[n].k = k + 1;
354:   s_u[n].c = 1;
355:   n++; /* Vy4 */
356:   s_u[n].i = i;
357:   s_u[n].j = j;
358:   s_u[n].k = k + 1;
359:   s_u[n].c = 2;
360:   n++; /* Vz4 */

362:   s_u[n].i = i;
363:   s_u[n].j = j + 1;
364:   s_u[n].k = k + 1;
365:   s_u[n].c = 0;
366:   n++;
367:   s_u[n].i = i;
368:   s_u[n].j = j + 1;
369:   s_u[n].k = k + 1;
370:   s_u[n].c = 1;
371:   n++;
372:   s_u[n].i = i;
373:   s_u[n].j = j + 1;
374:   s_u[n].k = k + 1;
375:   s_u[n].c = 2;
376:   n++;

378:   s_u[n].i = i + 1;
379:   s_u[n].j = j + 1;
380:   s_u[n].k = k + 1;
381:   s_u[n].c = 0;
382:   n++;
383:   s_u[n].i = i + 1;
384:   s_u[n].j = j + 1;
385:   s_u[n].k = k + 1;
386:   s_u[n].c = 1;
387:   n++;
388:   s_u[n].i = i + 1;
389:   s_u[n].j = j + 1;
390:   s_u[n].k = k + 1;
391:   s_u[n].c = 2;
392:   n++;

394:   s_u[n].i = i + 1;
395:   s_u[n].j = j;
396:   s_u[n].k = k + 1;
397:   s_u[n].c = 0;
398:   n++;
399:   s_u[n].i = i + 1;
400:   s_u[n].j = j;
401:   s_u[n].k = k + 1;
402:   s_u[n].c = 1;
403:   n++;
404:   s_u[n].i = i + 1;
405:   s_u[n].j = j;
406:   s_u[n].k = k + 1;
407:   s_u[n].c = 2;
408:   n++;

410:   /* pressure */
411:   n = 0;

413:   s_p[n].i = i;
414:   s_p[n].j = j;
415:   s_p[n].k = k;
416:   s_p[n].c = 3;
417:   n++; /* P0 */
418:   s_p[n].i = i;
419:   s_p[n].j = j + 1;
420:   s_p[n].k = k;
421:   s_p[n].c = 3;
422:   n++;
423:   s_p[n].i = i + 1;
424:   s_p[n].j = j + 1;
425:   s_p[n].k = k;
426:   s_p[n].c = 3;
427:   n++;
428:   s_p[n].i = i + 1;
429:   s_p[n].j = j;
430:   s_p[n].k = k;
431:   s_p[n].c = 3;
432:   n++;

434:   s_p[n].i = i;
435:   s_p[n].j = j;
436:   s_p[n].k = k + 1;
437:   s_p[n].c = 3;
438:   n++; /* P0 */
439:   s_p[n].i = i;
440:   s_p[n].j = j + 1;
441:   s_p[n].k = k + 1;
442:   s_p[n].c = 3;
443:   n++;
444:   s_p[n].i = i + 1;
445:   s_p[n].j = j + 1;
446:   s_p[n].k = k + 1;
447:   s_p[n].c = 3;
448:   n++;
449:   s_p[n].i = i + 1;
450:   s_p[n].j = j;
451:   s_p[n].k = k + 1;
452:   s_p[n].c = 3;
453:   n++;
454:   return 0;
455: }

457: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords, PetscInt i, PetscInt j, PetscInt k, PetscScalar el_coord[])
458: {
460:   /* get coords for the element */
461:   el_coord[0] = coords[k][j][i].x;
462:   el_coord[1] = coords[k][j][i].y;
463:   el_coord[2] = coords[k][j][i].z;

465:   el_coord[3] = coords[k][j + 1][i].x;
466:   el_coord[4] = coords[k][j + 1][i].y;
467:   el_coord[5] = coords[k][j + 1][i].z;

469:   el_coord[6] = coords[k][j + 1][i + 1].x;
470:   el_coord[7] = coords[k][j + 1][i + 1].y;
471:   el_coord[8] = coords[k][j + 1][i + 1].z;

473:   el_coord[9]  = coords[k][j][i + 1].x;
474:   el_coord[10] = coords[k][j][i + 1].y;
475:   el_coord[11] = coords[k][j][i + 1].z;

477:   el_coord[12] = coords[k + 1][j][i].x;
478:   el_coord[13] = coords[k + 1][j][i].y;
479:   el_coord[14] = coords[k + 1][j][i].z;

481:   el_coord[15] = coords[k + 1][j + 1][i].x;
482:   el_coord[16] = coords[k + 1][j + 1][i].y;
483:   el_coord[17] = coords[k + 1][j + 1][i].z;

485:   el_coord[18] = coords[k + 1][j + 1][i + 1].x;
486:   el_coord[19] = coords[k + 1][j + 1][i + 1].y;
487:   el_coord[20] = coords[k + 1][j + 1][i + 1].z;

489:   el_coord[21] = coords[k + 1][j][i + 1].x;
490:   el_coord[22] = coords[k + 1][j][i + 1].y;
491:   el_coord[23] = coords[k + 1][j][i + 1].z;
492:   return 0;
493: }

495: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field, PetscInt i, PetscInt j, PetscInt k, StokesDOF nodal_fields[])
496: {
498:   /* get the nodal fields for u */
499:   nodal_fields[0].u_dof = field[k][j][i].u_dof;
500:   nodal_fields[0].v_dof = field[k][j][i].v_dof;
501:   nodal_fields[0].w_dof = field[k][j][i].w_dof;

503:   nodal_fields[1].u_dof = field[k][j + 1][i].u_dof;
504:   nodal_fields[1].v_dof = field[k][j + 1][i].v_dof;
505:   nodal_fields[1].w_dof = field[k][j + 1][i].w_dof;

507:   nodal_fields[2].u_dof = field[k][j + 1][i + 1].u_dof;
508:   nodal_fields[2].v_dof = field[k][j + 1][i + 1].v_dof;
509:   nodal_fields[2].w_dof = field[k][j + 1][i + 1].w_dof;

511:   nodal_fields[3].u_dof = field[k][j][i + 1].u_dof;
512:   nodal_fields[3].v_dof = field[k][j][i + 1].v_dof;
513:   nodal_fields[3].w_dof = field[k][j][i + 1].w_dof;

515:   nodal_fields[4].u_dof = field[k + 1][j][i].u_dof;
516:   nodal_fields[4].v_dof = field[k + 1][j][i].v_dof;
517:   nodal_fields[4].w_dof = field[k + 1][j][i].w_dof;

519:   nodal_fields[5].u_dof = field[k + 1][j + 1][i].u_dof;
520:   nodal_fields[5].v_dof = field[k + 1][j + 1][i].v_dof;
521:   nodal_fields[5].w_dof = field[k + 1][j + 1][i].w_dof;

523:   nodal_fields[6].u_dof = field[k + 1][j + 1][i + 1].u_dof;
524:   nodal_fields[6].v_dof = field[k + 1][j + 1][i + 1].v_dof;
525:   nodal_fields[6].w_dof = field[k + 1][j + 1][i + 1].w_dof;

527:   nodal_fields[7].u_dof = field[k + 1][j][i + 1].u_dof;
528:   nodal_fields[7].v_dof = field[k + 1][j][i + 1].v_dof;
529:   nodal_fields[7].w_dof = field[k + 1][j][i + 1].w_dof;

531:   /* get the nodal fields for p */
532:   nodal_fields[0].p_dof = field[k][j][i].p_dof;
533:   nodal_fields[1].p_dof = field[k][j + 1][i].p_dof;
534:   nodal_fields[2].p_dof = field[k][j + 1][i + 1].p_dof;
535:   nodal_fields[3].p_dof = field[k][j][i + 1].p_dof;

537:   nodal_fields[4].p_dof = field[k + 1][j][i].p_dof;
538:   nodal_fields[5].p_dof = field[k + 1][j + 1][i].p_dof;
539:   nodal_fields[6].p_dof = field[k + 1][j + 1][i + 1].p_dof;
540:   nodal_fields[7].p_dof = field[k + 1][j][i + 1].p_dof;
541:   return 0;
542: }

544: 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)
545: {
546:   PetscInt              ij;
547:   PETSC_UNUSED PetscInt r, c, nr, nc;

549:   nr = w_NPE * w_dof;
550:   nc = u_NPE * u_dof;

552:   r = w_dof * wi + wd;
553:   c = u_dof * ui + ud;

555:   ij = r * nc + c;

557:   return ij;
558: }

560: static PetscErrorCode DMDASetValuesLocalStencil3D_ADD_VALUES(StokesDOF ***fields_F, MatStencil u_eqn[], MatStencil p_eqn[], PetscScalar Fe_u[], PetscScalar Fe_p[])
561: {
562:   PetscInt n, II, J, K;

565:   for (n = 0; n < NODES_PER_EL; n++) {
566:     II = u_eqn[NSD * n].i;
567:     J  = u_eqn[NSD * n].j;
568:     K  = u_eqn[NSD * n].k;

570:     fields_F[K][J][II].u_dof = fields_F[K][J][II].u_dof + Fe_u[NSD * n];

572:     II = u_eqn[NSD * n + 1].i;
573:     J  = u_eqn[NSD * n + 1].j;
574:     K  = u_eqn[NSD * n + 1].k;

576:     fields_F[K][J][II].v_dof = fields_F[K][J][II].v_dof + Fe_u[NSD * n + 1];

578:     II                       = u_eqn[NSD * n + 2].i;
579:     J                        = u_eqn[NSD * n + 2].j;
580:     K                        = u_eqn[NSD * n + 2].k;
581:     fields_F[K][J][II].w_dof = fields_F[K][J][II].w_dof + Fe_u[NSD * n + 2];

583:     II = p_eqn[n].i;
584:     J  = p_eqn[n].j;
585:     K  = p_eqn[n].k;

587:     fields_F[K][J][II].p_dof = fields_F[K][J][II].p_dof + Fe_p[n];
588:   }
589:   return 0;
590: }

592: static void FormStressOperatorQ13D(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
593: {
594:   PetscInt       ngp;
595:   PetscScalar    gp_xi[GAUSS_POINTS][NSD];
596:   PetscScalar    gp_weight[GAUSS_POINTS];
597:   PetscInt       p, i, j, k;
598:   PetscScalar    GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
599:   PetscScalar    J_p, tildeD[6];
600:   PetscScalar    B[6][U_DOFS * NODES_PER_EL];
601:   const PetscInt nvdof = U_DOFS * NODES_PER_EL;

603:   /* define quadrature rule */
604:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

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

611:     for (i = 0; i < NODES_PER_EL; i++) {
612:       PetscScalar d_dx_i = GNx_p[0][i];
613:       PetscScalar d_dy_i = GNx_p[1][i];
614:       PetscScalar d_dz_i = GNx_p[2][i];

616:       B[0][3 * i]     = d_dx_i;
617:       B[0][3 * i + 1] = 0.0;
618:       B[0][3 * i + 2] = 0.0;
619:       B[1][3 * i]     = 0.0;
620:       B[1][3 * i + 1] = d_dy_i;
621:       B[1][3 * i + 2] = 0.0;
622:       B[2][3 * i]     = 0.0;
623:       B[2][3 * i + 1] = 0.0;
624:       B[2][3 * i + 2] = d_dz_i;

626:       B[3][3 * i]     = d_dy_i;
627:       B[3][3 * i + 1] = d_dx_i;
628:       B[3][3 * i + 2] = 0.0; /* e_xy */
629:       B[4][3 * i]     = d_dz_i;
630:       B[4][3 * i + 1] = 0.0;
631:       B[4][3 * i + 2] = d_dx_i; /* e_xz */
632:       B[5][3 * i]     = 0.0;
633:       B[5][3 * i + 1] = d_dz_i;
634:       B[5][3 * i + 2] = d_dy_i; /* e_yz */
635:     }

637:     tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
638:     tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
639:     tildeD[2] = 2.0 * gp_weight[p] * J_p * eta[p];

641:     tildeD[3] = gp_weight[p] * J_p * eta[p];
642:     tildeD[4] = gp_weight[p] * J_p * eta[p];
643:     tildeD[5] = gp_weight[p] * J_p * eta[p];

645:     /* form Bt tildeD B */
646:     /*
647:      Ke_ij = Bt_ik . D_kl . B_lj
648:      = B_ki . D_kl . B_lj
649:      = B_ki . D_kk . B_kj
650:      */
651:     for (i = 0; i < nvdof; i++) {
652:       for (j = i; j < nvdof; j++) {
653:         for (k = 0; k < 6; k++) Ke[i * nvdof + j] += B[k][i] * tildeD[k] * B[k][j];
654:       }
655:     }
656:   }
657:   /* fill lower triangular part */
658: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
659:   for (i = 0; i < nvdof; i++) {
660:     for (j = i; j < nvdof; j++) Ke[j * nvdof + i] = Ke[i * nvdof + j];
661:   }
662: #endif
663: }

665: static void FormGradientOperatorQ13D(PetscScalar Ke[], PetscScalar coords[])
666: {
667:   PetscInt    ngp;
668:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
669:   PetscScalar gp_weight[GAUSS_POINTS];
670:   PetscInt    p, i, j, di;
671:   PetscScalar Ni_p[NODES_PER_EL];
672:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
673:   PetscScalar J_p, fac;

675:   /* define quadrature rule */
676:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

678:   /* evaluate integral */
679:   for (p = 0; p < ngp; p++) {
680:     ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
681:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
682:     ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
683:     fac = gp_weight[p] * J_p;

685:     for (i = 0; i < NODES_PER_EL; i++) {     /* u nodes */
686:       for (di = 0; di < NSD; di++) {         /* u dofs */
687:         for (j = 0; j < NODES_PER_EL; j++) { /* p nodes, p dofs = 1 (ie no loop) */
688:           PetscInt IJ;
689:           IJ = ASS_MAP_wIwDI_uJuDJ(i, di, NODES_PER_EL, 3, j, 0, NODES_PER_EL, 1);

691:           Ke[IJ] = Ke[IJ] - GNx_p[di][i] * Ni_p[j] * fac;
692:         }
693:       }
694:     }
695:   }
696: }

698: static void FormDivergenceOperatorQ13D(PetscScalar De[], PetscScalar coords[])
699: {
700:   PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
701:   PetscInt    i, j;
702:   PetscInt    nr_g, nc_g;

704:   PetscMemzero(Ge, sizeof(Ge));
705:   FormGradientOperatorQ13D(Ge, coords);

707:   nr_g = U_DOFS * NODES_PER_EL;
708:   nc_g = P_DOFS * NODES_PER_EL;

710:   for (i = 0; i < nr_g; i++) {
711:     for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
712:   }
713: }

715: static void FormStabilisationOperatorQ13D(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
716: {
717:   PetscInt    ngp;
718:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
719:   PetscScalar gp_weight[GAUSS_POINTS];
720:   PetscInt    p, i, j;
721:   PetscScalar Ni_p[NODES_PER_EL];
722:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
723:   PetscScalar J_p, fac, eta_avg;

725:   /* define quadrature rule */
726:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

728:   /* evaluate integral */
729:   for (p = 0; p < ngp; p++) {
730:     ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
731:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
732:     ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
733:     fac = gp_weight[p] * J_p;
734:     /*
735:      for (i = 0; i < NODES_PER_EL; i++) {
736:      for (j = i; j < NODES_PER_EL; j++) {
737:      Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
738:      }
739:      }
740:      */

742:     for (i = 0; i < NODES_PER_EL; i++) {
743:       for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * (Ni_p[i] * Ni_p[j] - 0.015625);
744:     }
745:   }

747:   /* scale */
748:   eta_avg = 0.0;
749:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
750:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
751:   fac     = 1.0 / eta_avg;

753:   /*
754:    for (i = 0; i < NODES_PER_EL; i++) {
755:    for (j = i; j < NODES_PER_EL; j++) {
756:    Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
757:    #if defined(ASSEMBLE_LOWER_TRIANGULAR)
758:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
759:    #endif
760:    }
761:    }
762:    */

764:   for (i = 0; i < NODES_PER_EL; i++) {
765:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
766:   }
767: }

769: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
770: {
771:   PetscInt    ngp;
772:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
773:   PetscScalar gp_weight[GAUSS_POINTS];
774:   PetscInt    p, i, j;
775:   PetscScalar Ni_p[NODES_PER_EL];
776:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
777:   PetscScalar J_p, fac, eta_avg;

779:   /* define quadrature rule */
780:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

782:   /* evaluate integral */
783:   for (p = 0; p < ngp; p++) {
784:     ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
785:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
786:     ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
787:     fac = gp_weight[p] * J_p;

789:     /*
790:      for (i = 0; i < NODES_PER_EL; i++) {
791:      for (j = i; j < NODES_PER_EL; j++) {
792:      Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
793:      }
794:      }
795:      */

797:     for (i = 0; i < NODES_PER_EL; i++) {
798:       for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = Ke[NODES_PER_EL * i + j] - fac * Ni_p[i] * Ni_p[j];
799:     }
800:   }

802:   /* scale */
803:   eta_avg = 0.0;
804:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
805:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
806:   fac     = 1.0 / eta_avg;
807:   /*
808:    for (i = 0; i < NODES_PER_EL; i++) {
809:    for (j = i; j < NODES_PER_EL; j++) {
810:    Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
811:    #if defined(ASSEMBLE_LOWER_TRIANGULAR)
812:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
813:    #endif
814:    }
815:    }
816:    */

818:   for (i = 0; i < NODES_PER_EL; i++) {
819:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
820:   }
821: }

823: static void FormMomentumRhsQ13D(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[], PetscScalar fz[])
824: {
825:   PetscInt    ngp;
826:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
827:   PetscScalar gp_weight[GAUSS_POINTS];
828:   PetscInt    p, i;
829:   PetscScalar Ni_p[NODES_PER_EL];
830:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
831:   PetscScalar J_p, fac;

833:   /* define quadrature rule */
834:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

836:   /* evaluate integral */
837:   for (p = 0; p < ngp; p++) {
838:     ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
839:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
840:     ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
841:     fac = gp_weight[p] * J_p;

843:     for (i = 0; i < NODES_PER_EL; i++) {
844:       Fe[NSD * i] -= fac * Ni_p[i] * fx[p];
845:       Fe[NSD * i + 1] -= fac * Ni_p[i] * fy[p];
846:       Fe[NSD * i + 2] -= fac * Ni_p[i] * fz[p];
847:     }
848:   }
849: }

851: static void FormContinuityRhsQ13D(PetscScalar Fe[], PetscScalar coords[], PetscScalar hc[])
852: {
853:   PetscInt    ngp;
854:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
855:   PetscScalar gp_weight[GAUSS_POINTS];
856:   PetscInt    p, i;
857:   PetscScalar Ni_p[NODES_PER_EL];
858:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
859:   PetscScalar J_p, fac;

861:   /* define quadrature rule */
862:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

864:   /* evaluate integral */
865:   for (p = 0; p < ngp; p++) {
866:     ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
867:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
868:     ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
869:     fac = gp_weight[p] * J_p;

871:     for (i = 0; i < NODES_PER_EL; i++) Fe[i] -= fac * Ni_p[i] * hc[p];
872:   }
873: }

875: #define _ZERO_ROWCOL_i(A, i) \
876:   { \
877:     PetscInt    KK; \
878:     PetscScalar tmp = A[24 * (i) + (i)]; \
879:     for (KK = 0; KK < 24; KK++) A[24 * (i) + KK] = 0.0; \
880:     for (KK = 0; KK < 24; KK++) A[24 * KK + (i)] = 0.0; \
881:     A[24 * (i) + (i)] = tmp; \
882:   }

884: #define _ZERO_ROW_i(A, i) \
885:   { \
886:     PetscInt KK; \
887:     for (KK = 0; KK < 8; KK++) A[8 * (i) + KK] = 0.0; \
888:   }

890: #define _ZERO_COL_i(A, i) \
891:   { \
892:     PetscInt KK; \
893:     for (KK = 0; KK < 8; KK++) A[24 * KK + (i)] = 0.0; \
894:   }

896: static PetscErrorCode AssembleA_Stokes(Mat A, DM stokes_da, CellProperties cell_properties)
897: {
898:   DM                      cda;
899:   Vec                     coords;
900:   DMDACoor3d           ***_coords;
901:   MatStencil              u_eqn[NODES_PER_EL * U_DOFS];
902:   MatStencil              p_eqn[NODES_PER_EL * P_DOFS];
903:   PetscInt                sex, sey, sez, mx, my, mz;
904:   PetscInt                ei, ej, ek;
905:   PetscScalar             Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
906:   PetscScalar             Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
907:   PetscScalar             De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
908:   PetscScalar             Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
909:   PetscScalar             el_coords[NODES_PER_EL * NSD];
910:   GaussPointCoefficients *props;
911:   PetscScalar            *prop_eta;
912:   PetscInt                n, M, N, P;

915:   DMDAGetInfo(stokes_da, 0, &M, &N, &P, 0, 0, 0, 0, 0, 0, 0, 0, 0);
916:   /* setup for coords */
917:   DMGetCoordinateDM(stokes_da, &cda);
918:   DMGetCoordinatesLocal(stokes_da, &coords);
919:   DMDAVecGetArray(cda, coords, &_coords);
920:   DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez);
921:   DMDAGetElementsSizes(stokes_da, &mx, &my, &mz);
922:   for (ek = sez; ek < sez + mz; ek++) {
923:     for (ej = sey; ej < sey + my; ej++) {
924:       for (ei = sex; ei < sex + mx; ei++) {
925:         /* get coords for the element */
926:         GetElementCoords3D(_coords, ei, ej, ek, el_coords);
927:         /* get cell properties */
928:         CellPropertiesGetCell(cell_properties, ei, ej, ek, &props);
929:         /* get coefficients for the element */
930:         prop_eta = props->eta;

932:         /* initialise element stiffness matrix */
933:         PetscMemzero(Ae, sizeof(Ae));
934:         PetscMemzero(Ge, sizeof(Ge));
935:         PetscMemzero(De, sizeof(De));
936:         PetscMemzero(Ce, sizeof(Ce));

938:         /* form element stiffness matrix */
939:         FormStressOperatorQ13D(Ae, el_coords, prop_eta);
940:         FormGradientOperatorQ13D(Ge, el_coords);
941:         /*#if defined(ASSEMBLE_LOWER_TRIANGULAR)*/
942:         FormDivergenceOperatorQ13D(De, el_coords);
943:         /*#endif*/
944:         FormStabilisationOperatorQ13D(Ce, el_coords, prop_eta);

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

949:         for (n = 0; n < NODES_PER_EL; n++) {
950:           if ((u_eqn[3 * n].i == 0) || (u_eqn[3 * n].i == M - 1)) {
951:             _ZERO_ROWCOL_i(Ae, 3 * n);
952:             _ZERO_ROW_i(Ge, 3 * n);
953:             _ZERO_COL_i(De, 3 * n);
954:           }

956:           if ((u_eqn[3 * n + 1].j == 0) || (u_eqn[3 * n + 1].j == N - 1)) {
957:             _ZERO_ROWCOL_i(Ae, 3 * n + 1);
958:             _ZERO_ROW_i(Ge, 3 * n + 1);
959:             _ZERO_COL_i(De, 3 * n + 1);
960:           }

962:           if ((u_eqn[3 * n + 2].k == 0) || (u_eqn[3 * n + 2].k == P - 1)) {
963:             _ZERO_ROWCOL_i(Ae, 3 * n + 2);
964:             _ZERO_ROW_i(Ge, 3 * n + 2);
965:             _ZERO_COL_i(De, 3 * n + 2);
966:           }
967:         }
968:         MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES);
969:         MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES);
970:         MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES);
971:         MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES);
972:       }
973:     }
974:   }
975:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
976:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

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

980:   return 0;
981: }

983: static PetscErrorCode AssembleA_PCStokes(Mat A, DM stokes_da, CellProperties cell_properties)
984: {
985:   DM                      cda;
986:   Vec                     coords;
987:   DMDACoor3d           ***_coords;
988:   MatStencil              u_eqn[NODES_PER_EL * U_DOFS];
989:   MatStencil              p_eqn[NODES_PER_EL * P_DOFS];
990:   PetscInt                sex, sey, sez, mx, my, mz;
991:   PetscInt                ei, ej, ek;
992:   PetscScalar             Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
993:   PetscScalar             Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
994:   PetscScalar             De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
995:   PetscScalar             Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
996:   PetscScalar             el_coords[NODES_PER_EL * NSD];
997:   GaussPointCoefficients *props;
998:   PetscScalar            *prop_eta;
999:   PetscInt                n, M, N, P;

1002:   DMDAGetInfo(stokes_da, 0, &M, &N, &P, 0, 0, 0, 0, 0, 0, 0, 0, 0);
1003:   /* setup for coords */
1004:   DMGetCoordinateDM(stokes_da, &cda);
1005:   DMGetCoordinatesLocal(stokes_da, &coords);
1006:   DMDAVecGetArray(cda, coords, &_coords);

1008:   DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez);
1009:   DMDAGetElementsSizes(stokes_da, &mx, &my, &mz);
1010:   for (ek = sez; ek < sez + mz; ek++) {
1011:     for (ej = sey; ej < sey + my; ej++) {
1012:       for (ei = sex; ei < sex + mx; ei++) {
1013:         /* get coords for the element */
1014:         GetElementCoords3D(_coords, ei, ej, ek, el_coords);
1015:         /* get cell properties */
1016:         CellPropertiesGetCell(cell_properties, ei, ej, ek, &props);
1017:         /* get coefficients for the element */
1018:         prop_eta = props->eta;

1020:         /* initialise element stiffness matrix */
1021:         PetscMemzero(Ae, sizeof(Ae));
1022:         PetscMemzero(Ge, sizeof(Ge));
1023:         PetscMemzero(De, sizeof(De));
1024:         PetscMemzero(Ce, sizeof(Ce));

1026:         /* form element stiffness matrix */
1027:         FormStressOperatorQ13D(Ae, el_coords, prop_eta);
1028:         FormGradientOperatorQ13D(Ge, el_coords);
1029:         /* FormDivergenceOperatorQ13D(De,el_coords); */
1030:         FormScaledMassMatrixOperatorQ13D(Ce, el_coords, prop_eta);

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

1035:         for (n = 0; n < NODES_PER_EL; n++) {
1036:           if ((u_eqn[3 * n].i == 0) || (u_eqn[3 * n].i == M - 1)) {
1037:             _ZERO_ROWCOL_i(Ae, 3 * n);
1038:             _ZERO_ROW_i(Ge, 3 * n);
1039:             _ZERO_COL_i(De, 3 * n);
1040:           }

1042:           if ((u_eqn[3 * n + 1].j == 0) || (u_eqn[3 * n + 1].j == N - 1)) {
1043:             _ZERO_ROWCOL_i(Ae, 3 * n + 1);
1044:             _ZERO_ROW_i(Ge, 3 * n + 1);
1045:             _ZERO_COL_i(De, 3 * n + 1);
1046:           }

1048:           if ((u_eqn[3 * n + 2].k == 0) || (u_eqn[3 * n + 2].k == P - 1)) {
1049:             _ZERO_ROWCOL_i(Ae, 3 * n + 2);
1050:             _ZERO_ROW_i(Ge, 3 * n + 2);
1051:             _ZERO_COL_i(De, 3 * n + 2);
1052:           }
1053:         }
1054:         MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES);
1055:         MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES);
1056:         /*MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);*/
1057:         MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES);
1058:       }
1059:     }
1060:   }
1061:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1062:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

1064:   DMDAVecRestoreArray(cda, coords, &_coords);
1065:   return 0;
1066: }

1068: static PetscErrorCode AssembleF_Stokes(Vec F, DM stokes_da, CellProperties cell_properties)
1069: {
1070:   DM                      cda;
1071:   Vec                     coords;
1072:   DMDACoor3d           ***_coords;
1073:   MatStencil              u_eqn[NODES_PER_EL * U_DOFS];
1074:   MatStencil              p_eqn[NODES_PER_EL * P_DOFS];
1075:   PetscInt                sex, sey, sez, mx, my, mz;
1076:   PetscInt                ei, ej, ek;
1077:   PetscScalar             Fe[NODES_PER_EL * U_DOFS];
1078:   PetscScalar             He[NODES_PER_EL * P_DOFS];
1079:   PetscScalar             el_coords[NODES_PER_EL * NSD];
1080:   GaussPointCoefficients *props;
1081:   PetscScalar            *prop_fx, *prop_fy, *prop_fz, *prop_hc;
1082:   Vec                     local_F;
1083:   StokesDOF            ***ff;
1084:   PetscInt                n, M, N, P;

1087:   DMDAGetInfo(stokes_da, 0, &M, &N, &P, 0, 0, 0, 0, 0, 0, 0, 0, 0);
1088:   /* setup for coords */
1089:   DMGetCoordinateDM(stokes_da, &cda);
1090:   DMGetCoordinatesLocal(stokes_da, &coords);
1091:   DMDAVecGetArray(cda, coords, &_coords);

1093:   /* get access to the vector */
1094:   DMGetLocalVector(stokes_da, &local_F);
1095:   VecZeroEntries(local_F);
1096:   DMDAVecGetArray(stokes_da, local_F, &ff);
1097:   DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez);
1098:   DMDAGetElementsSizes(stokes_da, &mx, &my, &mz);
1099:   for (ek = sez; ek < sez + mz; ek++) {
1100:     for (ej = sey; ej < sey + my; ej++) {
1101:       for (ei = sex; ei < sex + mx; ei++) {
1102:         /* get coords for the element */
1103:         GetElementCoords3D(_coords, ei, ej, ek, el_coords);
1104:         /* get cell properties */
1105:         CellPropertiesGetCell(cell_properties, ei, ej, ek, &props);
1106:         /* get coefficients for the element */
1107:         prop_fx = props->fx;
1108:         prop_fy = props->fy;
1109:         prop_fz = props->fz;
1110:         prop_hc = props->hc;

1112:         /* initialise element stiffness matrix */
1113:         PetscMemzero(Fe, sizeof(Fe));
1114:         PetscMemzero(He, sizeof(He));

1116:         /* form element stiffness matrix */
1117:         FormMomentumRhsQ13D(Fe, el_coords, prop_fx, prop_fy, prop_fz);
1118:         FormContinuityRhsQ13D(He, el_coords, prop_hc);

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

1123:         for (n = 0; n < NODES_PER_EL; n++) {
1124:           if ((u_eqn[3 * n].i == 0) || (u_eqn[3 * n].i == M - 1)) Fe[3 * n] = 0.0;

1126:           if ((u_eqn[3 * n + 1].j == 0) || (u_eqn[3 * n + 1].j == N - 1)) Fe[3 * n + 1] = 0.0;

1128:           if ((u_eqn[3 * n + 2].k == 0) || (u_eqn[3 * n + 2].k == P - 1)) Fe[3 * n + 2] = 0.0;
1129:         }

1131:         DMDASetValuesLocalStencil3D_ADD_VALUES(ff, u_eqn, p_eqn, Fe, He);
1132:       }
1133:     }
1134:   }
1135:   DMDAVecRestoreArray(stokes_da, local_F, &ff);
1136:   DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F);
1137:   DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F);
1138:   DMRestoreLocalVector(stokes_da, &local_F);

1140:   DMDAVecRestoreArray(cda, coords, &_coords);
1141:   return 0;
1142: }

1144: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta, PetscReal *MX, PetscReal *MY, PetscReal *MZ)
1145: {
1146:   *theta = 0.0;
1147:   *MX    = 2.0 * PETSC_PI;
1148:   *MY    = 2.0 * PETSC_PI;
1149:   *MZ    = 2.0 * PETSC_PI;
1150: }
1151: static void evaluate_MS_FrankKamentski(PetscReal pos[], PetscReal v[], PetscReal *p, PetscReal *eta, PetscReal Fm[], PetscReal *Fc)
1152: {
1153:   PetscReal x, y, z;
1154:   PetscReal theta, MX, MY, MZ;

1156:   evaluate_MS_FrankKamentski_constants(&theta, &MX, &MY, &MZ);
1157:   x = pos[0];
1158:   y = pos[1];
1159:   z = pos[2];
1160:   if (v) {
1161:     /*
1162:      v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1163:      v[1] = z*cos(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1164:      v[2] = -(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z);
1165:      */
1166:     /*
1167:      v[0] = PetscSinReal(PETSC_PI*x);
1168:      v[1] = PetscSinReal(PETSC_PI*y);
1169:      v[2] = PetscSinReal(PETSC_PI*z);
1170:      */
1171:     v[0] = PetscPowRealInt(z, 3) * PetscExpReal(y) * PetscSinReal(PETSC_PI * x);
1172:     v[1] = z * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y);
1173:     v[2] = PetscPowRealInt(z, 2) * (PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 - PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) / 2) - PETSC_PI * PetscPowRealInt(z, 4) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) / 4;
1174:   }
1175:   if (p) *p = PetscPowRealInt(x, 2) + PetscPowRealInt(y, 2) + PetscPowRealInt(z, 2);
1176:   if (eta) {
1177:     /* eta = PetscExpReal(-theta*(1.0 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));*/
1178:     *eta = 1.0;
1179:   }
1180:   if (Fm) {
1181:     /*
1182:      Fm[0] = -2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(PETSC_PI*x) - 0.2*MZ*theta*(-1.5*PetscPowRealInt(x,2)*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PetscPowRealInt(z,2)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x))*PetscCosReal(MX*x)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) - 0.2*PETSC_PI*MX*theta*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscCosReal(MZ*z)*PetscExpReal(y)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) + 2.0*(3.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 3.0*PETSC_PI*PetscPowRealInt(x,2)*PetscCosReal(2.0*PETSC_PI*z))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 1.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) ;
1183:      Fm[1] = -2*y - 0.2*MX*theta*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x))*PetscCosReal(MZ*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) - 0.2*MZ*theta*(-1.5*PetscPowRealInt(y,2)*PetscSinReal(2.0*PETSC_PI*z) + 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y))*PetscCosReal(MX*x)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) + 2.0*(-2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.5*PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 2*PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(-z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) - 6.0*PETSC_PI*PetscPowRealInt(y,2)*PetscCosReal(2.0*PETSC_PI*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));
1184:      Fm[2] = -2*z + 8.0*PetscPowRealInt(PETSC_PI,2)*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(2.0*PETSC_PI*z) - 0.2*MX*theta*(-1.5*PetscPowRealInt(x,2)*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PetscPowRealInt(z,2)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x))*PetscCosReal(MZ*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) + 0.4*PETSC_PI*MZ*theta*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(MX*x)*PetscCosReal(2.0*PETSC_PI*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) + 2.0*(-3.0*x*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(-3.0*y*PetscSinReal(2.0*PETSC_PI*z) - 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.5*PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(-1.5*PetscPowRealInt(y,2)*PetscSinReal(2.0*PETSC_PI*z) + 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))  ;
1185:      */
1186:     /*
1187:      Fm[0]=-2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*x);
1188:      Fm[1]=-2*y - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*y);
1189:      Fm[2]=-2*z - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*z);
1190:      */
1191:     /*
1192:      Fm[0] = -2*x + PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 6.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 6.0*PETSC_PI*PetscPowRealInt(x,2)*PetscCosReal(2.0*PETSC_PI*z) - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 2.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x) ;
1193:      Fm[1] = -2*y - 6.0*PETSC_PI*PetscPowRealInt(y,2)*PetscCosReal(2.0*PETSC_PI*z) + 2.0*z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 6.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) - 4.0*PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);
1194:      Fm[2] = -2*z - 6.0*x*PetscSinReal(2.0*PETSC_PI*z) - 6.0*y*PetscSinReal(2.0*PETSC_PI*z) - PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 8.0*PetscPowRealInt(PETSC_PI,2)*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z) + 3.0*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y) ;
1195:      */
1196:     Fm[0] = -2 * x + 2 * z * (PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(PETSC_PI * y) * PetscExpReal(-y) * PetscSinReal(2.0 * PETSC_PI * x) - 1.0 * PETSC_PI * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) * PetscSinReal(2.0 * PETSC_PI * x)) + PetscPowRealInt(z, 3) * PetscExpReal(y) * PetscSinReal(PETSC_PI * x) + 6.0 * z * PetscExpReal(y) * PetscSinReal(PETSC_PI * x) - 1.0 * PetscPowRealInt(PETSC_PI, 2) * PetscPowRealInt(z, 3) * PetscExpReal(y) * PetscSinReal(PETSC_PI * x) + 2.0 * PETSC_PI * z * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) * PetscSinReal(2.0 * PETSC_PI * x) - 2.0 * z * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(PETSC_PI * y) * PetscExpReal(-y) * PetscSinReal(2.0 * PETSC_PI * x);
1197:     Fm[1] = -2 * y + 2 * z * (-PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 + PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 + PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y)) + 2.0 * z * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) - 6.0 * z * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) - 4.0 * PETSC_PI * z * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y);
1198:     Fm[2] = -2 * z + PetscPowRealInt(z, 2) * (-2.0 * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) + 2.0 * PetscPowRealInt(PETSC_PI, 3) * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y)) + PetscPowRealInt(z, 2) * (PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 - 3 * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 + PetscPowRealInt(PETSC_PI, 3) * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) / 2 - 3 * PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) / 2) + 1.0 * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) + 0.25 * PetscPowRealInt(PETSC_PI, 3) * PetscPowRealInt(z, 4) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) - 0.25 * PETSC_PI * PetscPowRealInt(z, 4) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) - 3.0 * PETSC_PI * PetscPowRealInt(z, 2) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) - 1.0 * PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y);
1199:   }
1200:   if (Fc) {
1201:     /* Fc = -2.0*PETSC_PI*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(2.0*PETSC_PI*z) - z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y) ;*/
1202:     /* Fc = PETSC_PI*PetscCosReal(PETSC_PI*x) + PETSC_PI*PetscCosReal(PETSC_PI*y) + PETSC_PI*PetscCosReal(PETSC_PI*z);*/
1203:     /* Fc = -2.0*PETSC_PI*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(2.0*PETSC_PI*z) - z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);*/
1204:     *Fc = 0.0;
1205:   }
1206: }

1208: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx, PetscInt my, PetscInt mz, DM *_da, Vec *_X)
1209: {
1210:   DM            da, cda;
1211:   Vec           X;
1212:   StokesDOF  ***_stokes;
1213:   Vec           coords;
1214:   DMDACoor3d ***_coords;
1215:   PetscInt      si, sj, sk, ei, ej, ek, i, j, k;

1218:   DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, mz + 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 4, 1, NULL, NULL, NULL, &da);
1219:   DMSetFromOptions(da);
1220:   DMSetUp(da);
1221:   DMDASetFieldName(da, 0, "anlytic_Vx");
1222:   DMDASetFieldName(da, 1, "anlytic_Vy");
1223:   DMDASetFieldName(da, 2, "anlytic_Vz");
1224:   DMDASetFieldName(da, 3, "analytic_P");

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

1228:   DMGetCoordinatesLocal(da, &coords);
1229:   DMGetCoordinateDM(da, &cda);
1230:   DMDAVecGetArray(cda, coords, &_coords);

1232:   DMCreateGlobalVector(da, &X);
1233:   DMDAVecGetArray(da, X, &_stokes);

1235:   DMDAGetCorners(da, &si, &sj, &sk, &ei, &ej, &ek);
1236:   for (k = sk; k < sk + ek; k++) {
1237:     for (j = sj; j < sj + ej; j++) {
1238:       for (i = si; i < si + ei; i++) {
1239:         PetscReal pos[NSD], pressure, vel[NSD];

1241:         pos[0] = PetscRealPart(_coords[k][j][i].x);
1242:         pos[1] = PetscRealPart(_coords[k][j][i].y);
1243:         pos[2] = PetscRealPart(_coords[k][j][i].z);

1245:         evaluate_MS_FrankKamentski(pos, vel, &pressure, NULL, NULL, NULL);

1247:         _stokes[k][j][i].u_dof = vel[0];
1248:         _stokes[k][j][i].v_dof = vel[1];
1249:         _stokes[k][j][i].w_dof = vel[2];
1250:         _stokes[k][j][i].p_dof = pressure;
1251:       }
1252:     }
1253:   }
1254:   DMDAVecRestoreArray(da, X, &_stokes);
1255:   DMDAVecRestoreArray(cda, coords, &_coords);

1257:   *_da = da;
1258:   *_X  = X;
1259:   return 0;
1260: }

1262: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da, Vec X, Vec X_analytic)
1263: {
1264:   DM            cda;
1265:   Vec           coords, X_analytic_local, X_local;
1266:   DMDACoor3d ***_coords;
1267:   PetscInt      sex, sey, sez, mx, my, mz;
1268:   PetscInt      ei, ej, ek;
1269:   PetscScalar   el_coords[NODES_PER_EL * NSD];
1270:   StokesDOF  ***stokes_analytic, ***stokes;
1271:   StokesDOF     stokes_analytic_e[NODES_PER_EL], stokes_e[NODES_PER_EL];

1273:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
1274:   PetscScalar Ni_p[NODES_PER_EL];
1275:   PetscInt    ngp;
1276:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
1277:   PetscScalar gp_weight[GAUSS_POINTS];
1278:   PetscInt    p, i;
1279:   PetscScalar J_p, fac;
1280:   PetscScalar h, p_e_L2, u_e_L2, u_e_H1, p_L2, u_L2, u_H1, tp_L2, tu_L2, tu_H1;
1281:   PetscScalar tint_p_ms, tint_p, int_p_ms, int_p;
1282:   PetscInt    M;
1283:   PetscReal   xymin[NSD], xymax[NSD];

1286:   /* define quadrature rule */
1287:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

1289:   /* setup for coords */
1290:   DMGetCoordinateDM(stokes_da, &cda);
1291:   DMGetCoordinatesLocal(stokes_da, &coords);
1292:   DMDAVecGetArray(cda, coords, &_coords);

1294:   /* setup for analytic */
1295:   DMCreateLocalVector(stokes_da, &X_analytic_local);
1296:   DMGlobalToLocalBegin(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local);
1297:   DMGlobalToLocalEnd(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local);
1298:   DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic);

1300:   /* setup for solution */
1301:   DMCreateLocalVector(stokes_da, &X_local);
1302:   DMGlobalToLocalBegin(stokes_da, X, INSERT_VALUES, X_local);
1303:   DMGlobalToLocalEnd(stokes_da, X, INSERT_VALUES, X_local);
1304:   DMDAVecGetArray(stokes_da, X_local, &stokes);

1306:   DMDAGetInfo(stokes_da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
1307:   DMGetBoundingBox(stokes_da, xymin, xymax);

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

1311:   tp_L2 = tu_L2 = tu_H1 = 0.0;
1312:   tint_p_ms = tint_p = 0.0;

1314:   DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez);
1315:   DMDAGetElementsSizes(stokes_da, &mx, &my, &mz);
1316:   for (ek = sez; ek < sez + mz; ek++) {
1317:     for (ej = sey; ej < sey + my; ej++) {
1318:       for (ei = sex; ei < sex + mx; ei++) {
1319:         /* get coords for the element */
1320:         GetElementCoords3D(_coords, ei, ej, ek, el_coords);
1321:         StokesDAGetNodalFields3D(stokes, ei, ej, ek, stokes_e);
1322:         StokesDAGetNodalFields3D(stokes_analytic, ei, ej, ek, stokes_analytic_e);

1324:         /* evaluate integral */
1325:         for (p = 0; p < ngp; p++) {
1326:           ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
1327:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
1328:           ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, el_coords, &J_p);
1329:           fac = gp_weight[p] * J_p;

1331:           for (i = 0; i < NODES_PER_EL; i++) {
1332:             tint_p_ms = tint_p_ms + fac * Ni_p[i] * stokes_analytic_e[i].p_dof;
1333:             tint_p    = tint_p + fac * Ni_p[i] * stokes_e[i].p_dof;
1334:           }
1335:         }
1336:       }
1337:     }
1338:   }
1339:   MPI_Allreduce(&tint_p_ms, &int_p_ms, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD);
1340:   MPI_Allreduce(&tint_p, &int_p, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD);

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

1344:   /* remove mine and add manufacture one */
1345:   DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic);
1346:   DMDAVecRestoreArray(stokes_da, X_local, &stokes);

1348:   {
1349:     PetscInt     k, L, dof;
1350:     PetscScalar *fields;
1351:     DMDAGetInfo(stokes_da, 0, 0, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);

1353:     VecGetLocalSize(X_local, &L);
1354:     VecGetArray(X_local, &fields);
1355:     for (k = 0; k < L / dof; k++) fields[dof * k + 3] = fields[dof * k + 3] - int_p + int_p_ms;
1356:     VecRestoreArray(X_local, &fields);

1358:     VecGetLocalSize(X, &L);
1359:     VecGetArray(X, &fields);
1360:     for (k = 0; k < L / dof; k++) fields[dof * k + 3] = fields[dof * k + 3] - int_p + int_p_ms;
1361:     VecRestoreArray(X, &fields);
1362:   }

1364:   DMDAVecGetArray(stokes_da, X_local, &stokes);
1365:   DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic);

1367:   for (ek = sez; ek < sez + mz; ek++) {
1368:     for (ej = sey; ej < sey + my; ej++) {
1369:       for (ei = sex; ei < sex + mx; ei++) {
1370:         /* get coords for the element */
1371:         GetElementCoords3D(_coords, ei, ej, ek, el_coords);
1372:         StokesDAGetNodalFields3D(stokes, ei, ej, ek, stokes_e);
1373:         StokesDAGetNodalFields3D(stokes_analytic, ei, ej, ek, stokes_analytic_e);

1375:         /* evaluate integral */
1376:         p_e_L2 = 0.0;
1377:         u_e_L2 = 0.0;
1378:         u_e_H1 = 0.0;
1379:         for (p = 0; p < ngp; p++) {
1380:           ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
1381:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
1382:           ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, el_coords, &J_p);
1383:           fac = gp_weight[p] * J_p;

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

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

1390:             u_error = stokes_e[i].u_dof - stokes_analytic_e[i].u_dof;
1391:             v_error = stokes_e[i].v_dof - stokes_analytic_e[i].v_dof;
1392:             w_error = stokes_e[i].w_dof - stokes_analytic_e[i].w_dof;
1393:             u_e_L2 += fac * Ni_p[i] * (u_error * u_error + v_error * v_error + w_error * w_error);

1395:             u_e_H1 = u_e_H1 + fac * (GNx_p[0][i] * u_error * GNx_p[0][i] * u_error                                                                                                   /* du/dx */
1396:                                      + GNx_p[1][i] * u_error * GNx_p[1][i] * u_error + GNx_p[2][i] * u_error * GNx_p[2][i] * u_error + GNx_p[0][i] * v_error * GNx_p[0][i] * v_error /* dv/dx */
1397:                                      + GNx_p[1][i] * v_error * GNx_p[1][i] * v_error + GNx_p[2][i] * v_error * GNx_p[2][i] * v_error + GNx_p[0][i] * w_error * GNx_p[0][i] * w_error /* dw/dx */
1398:                                      + GNx_p[1][i] * w_error * GNx_p[1][i] * w_error + GNx_p[2][i] * w_error * GNx_p[2][i] * w_error);
1399:           }
1400:         }

1402:         tp_L2 += p_e_L2;
1403:         tu_L2 += u_e_L2;
1404:         tu_H1 += u_e_H1;
1405:       }
1406:     }
1407:   }
1408:   MPI_Allreduce(&tp_L2, &p_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD);
1409:   MPI_Allreduce(&tu_L2, &u_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD);
1410:   MPI_Allreduce(&tu_H1, &u_H1, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD);
1411:   p_L2 = PetscSqrtScalar(p_L2);
1412:   u_L2 = PetscSqrtScalar(u_L2);
1413:   u_H1 = PetscSqrtScalar(u_H1);

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

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

1419:   DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic);
1420:   VecDestroy(&X_analytic_local);
1421:   DMDAVecRestoreArray(stokes_da, X_local, &stokes);
1422:   VecDestroy(&X_local);
1423:   return 0;
1424: }

1426: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da, Vec FIELD, const char file_prefix[])
1427: {
1428:   char         vtk_filename[PETSC_MAX_PATH_LEN];
1429:   PetscMPIInt  rank;
1430:   MPI_Comm     comm;
1431:   FILE        *vtk_fp     = NULL;
1432:   const char  *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
1433:   PetscInt     si, sj, sk, nx, ny, nz, i;
1434:   PetscInt     f, n_fields, N;
1435:   DM           cda;
1436:   Vec          coords;
1437:   Vec          l_FIELD;
1438:   PetscScalar *_L_FIELD;
1439:   PetscInt     memory_offset;
1440:   PetscScalar *buffer;


1444:   /* create file name */
1445:   PetscObjectGetComm((PetscObject)da, &comm);
1446:   MPI_Comm_rank(comm, &rank);
1447:   PetscSNPrintf(vtk_filename, sizeof(vtk_filename), "subdomain-%s-p%1.4d.vts", file_prefix, rank);

1449:   /* open file and write header */
1450:   vtk_fp = fopen(vtk_filename, "w");

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

1455:   /* coords */
1456:   DMDAGetGhostCorners(da, &si, &sj, &sk, &nx, &ny, &nz);
1457:   N = nx * ny * nz;

1459:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order);
1460:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "  <StructuredGrid WholeExtent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", si, si + nx - 1, sj, sj + ny - 1, sk, sk + nz - 1);
1461:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    <Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", si, si + nx - 1, sj, sj + ny - 1, sk, sk + nz - 1);

1463:   memory_offset = 0;

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

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

1469:   /* copy coordinates */
1470:   DMGetCoordinateDM(da, &cda);
1471:   DMGetCoordinatesLocal(da, &coords);
1472:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "        <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n", memory_offset);
1473:   memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar) * N * 3;

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

1477:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "      <PointData Scalars=\" ");
1478:   DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_fields, 0, 0, 0, 0, 0);
1479:   for (f = 0; f < n_fields; f++) {
1480:     const char *field_name;
1481:     DMDAGetFieldName(da, f, &field_name);
1482:     PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "%s ", field_name);
1483:   }
1484:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "\">\n");

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

1489:     DMDAGetFieldName(da, f, &field_name);
1490:     PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "        <DataArray type=\"Float64\" Name=\"%s\" format=\"appended\" offset=\"%" PetscInt_FMT "\"/>\n", field_name, memory_offset);
1491:     memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar) * N;
1492:   }

1494:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "      </PointData>\n");
1495:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    </Piece>\n");
1496:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "  </StructuredGrid>\n");

1498:   PetscMalloc1(N, &buffer);
1499:   DMGetLocalVector(da, &l_FIELD);
1500:   DMGlobalToLocalBegin(da, FIELD, INSERT_VALUES, l_FIELD);
1501:   DMGlobalToLocalEnd(da, FIELD, INSERT_VALUES, l_FIELD);
1502:   VecGetArray(l_FIELD, &_L_FIELD);

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

1507:   /* write coordinates */
1508:   {
1509:     int          length = sizeof(PetscScalar) * N * 3;
1510:     PetscScalar *allcoords;

1512:     fwrite(&length, sizeof(int), 1, vtk_fp);
1513:     VecGetArray(coords, &allcoords);
1514:     fwrite(allcoords, sizeof(PetscScalar), 3 * N, vtk_fp);
1515:     VecRestoreArray(coords, &allcoords);
1516:   }
1517:   /* write fields */
1518:   for (f = 0; f < n_fields; f++) {
1519:     int length = sizeof(PetscScalar) * N;
1520:     fwrite(&length, sizeof(int), 1, vtk_fp);
1521:     /* load */
1522:     for (i = 0; i < N; i++) buffer[i] = _L_FIELD[n_fields * i + f];

1524:     /* write */
1525:     fwrite(buffer, sizeof(PetscScalar), N, vtk_fp);
1526:   }
1527:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "\n  </AppendedData>\n");

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

1531:   PetscFree(buffer);
1532:   VecRestoreArray(l_FIELD, &_L_FIELD);
1533:   DMRestoreLocalVector(da, &l_FIELD);

1535:   if (vtk_fp) {
1536:     fclose(vtk_fp);
1537:     vtk_fp = NULL;
1538:   }

1540:   return 0;
1541: }

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

1553:   /* create file name */
1554:   PetscObjectGetComm((PetscObject)da, &comm);
1555:   MPI_Comm_size(comm, &size);
1556:   MPI_Comm_rank(comm, &rank);

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

1561:   /* generate start,end list */
1562:   PetscMalloc1(pM + 1, &olx);
1563:   PetscMalloc1(pN + 1, &oly);
1564:   PetscMalloc1(pP + 1, &olz);
1565:   sum = 0;
1566:   for (i = 0; i < pM; i++) {
1567:     olx[i] = sum;
1568:     sum    = sum + lx[i];
1569:   }
1570:   olx[pM] = sum;
1571:   sum     = 0;
1572:   for (i = 0; i < pN; i++) {
1573:     oly[i] = sum;
1574:     sum    = sum + ly[i];
1575:   }
1576:   oly[pN] = sum;
1577:   sum     = 0;
1578:   for (i = 0; i < pP; i++) {
1579:     olz[i] = sum;
1580:     sum    = sum + lz[i];
1581:   }
1582:   olz[pP] = sum;

1584:   PetscMalloc1(pM, &osx);
1585:   PetscMalloc1(pN, &osy);
1586:   PetscMalloc1(pP, &osz);
1587:   PetscMalloc1(pM, &oex);
1588:   PetscMalloc1(pN, &oey);
1589:   PetscMalloc1(pP, &oez);
1590:   for (i = 0; i < pM; i++) {
1591:     osx[i] = olx[i] - stencil;
1592:     oex[i] = olx[i] + lx[i] + stencil;
1593:     if (osx[i] < 0) osx[i] = 0;
1594:     if (oex[i] > M) oex[i] = M;
1595:   }

1597:   for (i = 0; i < pN; i++) {
1598:     osy[i] = oly[i] - stencil;
1599:     oey[i] = oly[i] + ly[i] + stencil;
1600:     if (osy[i] < 0) osy[i] = 0;
1601:     if (oey[i] > M) oey[i] = N;
1602:   }
1603:   for (i = 0; i < pP; i++) {
1604:     osz[i] = olz[i] - stencil;
1605:     oez[i] = olz[i] + lz[i] + stencil;
1606:     if (osz[i] < 0) osz[i] = 0;
1607:     if (oez[i] > P) oez[i] = P;
1608:   }

1610:   for (k = 0; k < pP; k++) {
1611:     for (j = 0; j < pN; j++) {
1612:       for (i = 0; i < pM; i++) {
1613:         char     name[PETSC_MAX_PATH_LEN];
1614:         PetscInt procid = i + j * pM + k * pM * pN; /* convert proc(i,j,k) to pid */
1615:         PetscSNPrintf(name, sizeof(name), "subdomain-%s-p%1.4" PetscInt_FMT ".vts", local_file_prefix, procid);
1616:         for (II = 0; II < indent_level; II++) PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "  ");

1618:         PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\"      Source=\"%s\"/>\n", osx[i], oex[i] - 1, osy[j], oey[j] - 1, osz[k], oez[k] - 1, name);
1619:       }
1620:     }
1621:   }
1622:   PetscFree(olx);
1623:   PetscFree(oly);
1624:   PetscFree(olz);
1625:   PetscFree(osx);
1626:   PetscFree(osy);
1627:   PetscFree(osz);
1628:   PetscFree(oex);
1629:   PetscFree(oey);
1630:   PetscFree(oez);
1631:   return 0;
1632: }

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

1645:   /* only rank-0 generates this file */
1646:   PetscObjectGetComm((PetscObject)da, &comm);
1647:   MPI_Comm_size(comm, &size);
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");

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

1660:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order);

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

1667:   /* DUMP THE CELL REFERENCES */
1668:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    <PCellData>\n");
1669:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    </PCellData>\n");

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

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

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

1686:   /* close the file */
1687:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "  </PStructuredGrid>\n");
1688:   PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "</VTKFile>\n");

1690:   if (vtk_fp) {
1691:     fclose(vtk_fp);
1692:     vtk_fp = NULL;
1693:   }
1694:   return 0;
1695: }

1697: PetscErrorCode DAView3DPVTS(DM da, Vec x, const char NAME[])
1698: {
1699:   char vts_filename[PETSC_MAX_PATH_LEN];
1700:   char pvts_filename[PETSC_MAX_PATH_LEN];

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

1706:   PetscSNPrintf(pvts_filename, sizeof(pvts_filename), "%s-mesh", NAME);
1707:   DAView_3DVTK_PStructuredGrid(da, pvts_filename, vts_filename);
1708:   return 0;
1709: }

1711: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp, PetscInt n, PetscReal rnorm, void *dummy)
1712: {
1713:   PetscReal norms[4];
1714:   Vec       Br, v, w;
1715:   Mat       A;

1718:   KSPGetOperators(ksp, &A, NULL);
1719:   MatCreateVecs(A, &w, &v);

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

1723:   VecStrideNorm(Br, 0, NORM_2, &norms[0]);
1724:   VecStrideNorm(Br, 1, NORM_2, &norms[1]);
1725:   VecStrideNorm(Br, 2, NORM_2, &norms[2]);
1726:   VecStrideNorm(Br, 3, NORM_2, &norms[3]);

1728:   VecDestroy(&v);
1729:   VecDestroy(&w);

1731:   PetscPrintf(PETSC_COMM_WORLD, "%3" PetscInt_FMT " KSP Component U,V,W,P residual norm [ %1.12e, %1.12e, %1.12e, %1.12e ]\n", n, (double)norms[0], (double)norms[1], (double)norms[2], (double)norms[3]);
1732:   return 0;
1733: }

1735: static PetscErrorCode PCMGSetupViaCoarsen(PC pc, DM da_fine)
1736: {
1737:   PetscInt              nlevels, k;
1738:   PETSC_UNUSED PetscInt finest;
1739:   DM                   *da_list, *daclist;
1740:   Mat                   R;

1743:   nlevels = 1;
1744:   PetscOptionsGetInt(NULL, NULL, "-levels", &nlevels, 0);

1746:   PetscMalloc1(nlevels, &da_list);
1747:   for (k = 0; k < nlevels; k++) da_list[k] = NULL;
1748:   PetscMalloc1(nlevels, &daclist);
1749:   for (k = 0; k < nlevels; k++) daclist[k] = NULL;

1751:   /* finest grid is nlevels - 1 */
1752:   finest     = nlevels - 1;
1753:   daclist[0] = da_fine;
1754:   PetscObjectReference((PetscObject)da_fine);
1755:   DMCoarsenHierarchy(da_fine, nlevels - 1, &daclist[1]);
1756:   for (k = 0; k < nlevels; k++) {
1757:     da_list[k] = daclist[nlevels - 1 - k];
1758:     DMDASetUniformCoordinates(da_list[k], 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
1759:   }

1761:   PCMGSetLevels(pc, nlevels, NULL);
1762:   PCMGSetType(pc, PC_MG_MULTIPLICATIVE);
1763:   PCMGSetGalerkin(pc, PC_MG_GALERKIN_PMAT);

1765:   for (k = 1; k < nlevels; k++) {
1766:     DMCreateInterpolation(da_list[k - 1], da_list[k], &R, NULL);
1767:     PCMGSetInterpolation(pc, k, R);
1768:     MatDestroy(&R);
1769:   }

1771:   /* tidy up */
1772:   for (k = 0; k < nlevels; k++) DMDestroy(&da_list[k]);
1773:   PetscFree(da_list);
1774:   PetscFree(daclist);
1775:   return 0;
1776: }

1778: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx, PetscInt my, PetscInt mz)
1779: {
1780:   DM             da_Stokes;
1781:   PetscInt       u_dof, p_dof, dof, stencil_width;
1782:   Mat            A, B;
1783:   DM             vel_cda;
1784:   Vec            vel_coords;
1785:   PetscInt       p;
1786:   Vec            f, X, X1;
1787:   DMDACoor3d  ***_vel_coords;
1788:   PetscInt       its;
1789:   KSP            ksp_S;
1790:   PetscInt       model_definition = 0;
1791:   PetscInt       ei, ej, ek, sex, sey, sez, Imx, Jmy, Kmz;
1792:   CellProperties cell_properties;
1793:   PetscBool      write_output = PETSC_FALSE, resolve = PETSC_FALSE;

1796:   PetscOptionsGetBool(NULL, NULL, "-resolve", &resolve, NULL);
1797:   /* Generate the da for velocity and pressure */
1798:   /* Num nodes in each direction is mx+1, my+1, mz+1 */
1799:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1800:   p_dof         = P_DOFS; /* p - pressure */
1801:   dof           = u_dof + p_dof;
1802:   stencil_width = 1;
1803:   DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, mz + 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, NULL, &da_Stokes);
1804:   DMSetMatType(da_Stokes, MATAIJ);
1805:   DMSetFromOptions(da_Stokes);
1806:   DMSetUp(da_Stokes);
1807:   DMDASetFieldName(da_Stokes, 0, "Vx");
1808:   DMDASetFieldName(da_Stokes, 1, "Vy");
1809:   DMDASetFieldName(da_Stokes, 2, "Vz");
1810:   DMDASetFieldName(da_Stokes, 3, "P");

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

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

1818:   /* interpolate the coordinates to quadrature points */
1819:   DMGetCoordinateDM(da_Stokes, &vel_cda);
1820:   DMGetCoordinatesLocal(da_Stokes, &vel_coords);
1821:   DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords);
1822:   DMDAGetElementsCorners(da_Stokes, &sex, &sey, &sez);
1823:   DMDAGetElementsSizes(da_Stokes, &Imx, &Jmy, &Kmz);
1824:   for (ek = sez; ek < sez + Kmz; ek++) {
1825:     for (ej = sey; ej < sey + Jmy; ej++) {
1826:       for (ei = sex; ei < sex + Imx; ei++) {
1827:         /* get coords for the element */
1828:         PetscInt                ngp;
1829:         PetscScalar             gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
1830:         PetscScalar             el_coords[NSD * NODES_PER_EL];
1831:         GaussPointCoefficients *cell;

1833:         CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell);
1834:         GetElementCoords3D(_vel_coords, ei, ej, ek, el_coords);
1835:         ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

1837:         for (p = 0; p < GAUSS_POINTS; p++) {
1838:           PetscScalar xi_p[NSD], Ni_p[NODES_PER_EL];
1839:           PetscScalar gp_x, gp_y, gp_z;
1840:           PetscInt    n;

1842:           xi_p[0] = gp_xi[p][0];
1843:           xi_p[1] = gp_xi[p][1];
1844:           xi_p[2] = gp_xi[p][2];
1845:           ShapeFunctionQ13D_Evaluate(xi_p, Ni_p);

1847:           gp_x = gp_y = gp_z = 0.0;
1848:           for (n = 0; n < NODES_PER_EL; n++) {
1849:             gp_x = gp_x + Ni_p[n] * el_coords[NSD * n];
1850:             gp_y = gp_y + Ni_p[n] * el_coords[NSD * n + 1];
1851:             gp_z = gp_z + Ni_p[n] * el_coords[NSD * n + 2];
1852:           }
1853:           cell->gp_coords[NSD * p]     = gp_x;
1854:           cell->gp_coords[NSD * p + 1] = gp_y;
1855:           cell->gp_coords[NSD * p + 2] = gp_z;
1856:         }
1857:       }
1858:     }
1859:   }

1861:   PetscOptionsGetInt(NULL, NULL, "-model", &model_definition, NULL);

1863:   switch (model_definition) {
1864:   case 0: /* isoviscous */
1865:     for (ek = sez; ek < sez + Kmz; ek++) {
1866:       for (ej = sey; ej < sey + Jmy; ej++) {
1867:         for (ei = sex; ei < sex + Imx; ei++) {
1868:           GaussPointCoefficients *cell;

1870:           CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell);
1871:           for (p = 0; p < GAUSS_POINTS; p++) {
1872:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD * p]);
1873:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD * p + 1]);
1874:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD * p + 2]);

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

1878:             cell->fx[p] = 0.0 * coord_x;
1879:             cell->fy[p] = -PetscSinReal(2.2 * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1880:             cell->fz[p] = 0.0 * coord_z;
1881:             cell->hc[p] = 0.0;
1882:           }
1883:         }
1884:       }
1885:     }
1886:     break;

1888:   case 1: /* manufactured */
1889:     for (ek = sez; ek < sez + Kmz; ek++) {
1890:       for (ej = sey; ej < sey + Jmy; ej++) {
1891:         for (ei = sex; ei < sex + Imx; ei++) {
1892:           PetscReal               eta, Fm[NSD], Fc, pos[NSD];
1893:           GaussPointCoefficients *cell;

1895:           CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell);
1896:           for (p = 0; p < GAUSS_POINTS; p++) {
1897:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD * p]);
1898:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD * p + 1]);
1899:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD * p + 2]);

1901:             pos[0] = coord_x;
1902:             pos[1] = coord_y;
1903:             pos[2] = coord_z;

1905:             evaluate_MS_FrankKamentski(pos, NULL, NULL, &eta, Fm, &Fc);
1906:             cell->eta[p] = eta;

1908:             cell->fx[p] = Fm[0];
1909:             cell->fy[p] = Fm[1];
1910:             cell->fz[p] = Fm[2];
1911:             cell->hc[p] = Fc;
1912:           }
1913:         }
1914:       }
1915:     }
1916:     break;

1918:   case 2: /* solcx */
1919:     for (ek = sez; ek < sez + Kmz; ek++) {
1920:       for (ej = sey; ej < sey + Jmy; ej++) {
1921:         for (ei = sex; ei < sex + Imx; ei++) {
1922:           GaussPointCoefficients *cell;

1924:           CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell);
1925:           for (p = 0; p < GAUSS_POINTS; p++) {
1926:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD * p]);
1927:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD * p + 1]);
1928:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD * p + 2]);

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

1932:             cell->fx[p] = 0.0;
1933:             cell->fy[p] = -PetscSinReal(3.0 * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1934:             cell->fz[p] = 0.0 * coord_z;
1935:             cell->hc[p] = 0.0;
1936:           }
1937:         }
1938:       }
1939:     }
1940:     break;

1942:   case 3: /* sinker */
1943:     for (ek = sez; ek < sez + Kmz; ek++) {
1944:       for (ej = sey; ej < sey + Jmy; ej++) {
1945:         for (ei = sex; ei < sex + Imx; ei++) {
1946:           GaussPointCoefficients *cell;

1948:           CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell);
1949:           for (p = 0; p < GAUSS_POINTS; p++) {
1950:             PetscReal xp = PetscRealPart(cell->gp_coords[NSD * p]);
1951:             PetscReal yp = PetscRealPart(cell->gp_coords[NSD * p + 1]);
1952:             PetscReal zp = PetscRealPart(cell->gp_coords[NSD * p + 2]);

1954:             cell->eta[p] = 1.0e-2;
1955:             cell->fx[p]  = 0.0;
1956:             cell->fy[p]  = 0.0;
1957:             cell->fz[p]  = 0.0;
1958:             cell->hc[p]  = 0.0;

1960:             if ((PetscAbs(xp - 0.5) < 0.2) && (PetscAbs(yp - 0.5) < 0.2) && (PetscAbs(zp - 0.5) < 0.2)) {
1961:               cell->eta[p] = 1.0;
1962:               cell->fz[p]  = 1.0;
1963:             }
1964:           }
1965:         }
1966:       }
1967:     }
1968:     break;

1970:   case 4: /* subdomain jumps */
1971:     for (ek = sez; ek < sez + Kmz; ek++) {
1972:       for (ej = sey; ej < sey + Jmy; ej++) {
1973:         for (ei = sex; ei < sex + Imx; ei++) {
1974:           PetscReal               opts_mag, opts_eta0;
1975:           PetscInt                px, py, pz;
1976:           PetscBool               jump;
1977:           PetscMPIInt             rr;
1978:           GaussPointCoefficients *cell;

1980:           opts_mag  = 1.0;
1981:           opts_eta0 = 1.e-2;

1983:           PetscOptionsGetReal(NULL, NULL, "-jump_eta0", &opts_eta0, NULL);
1984:           PetscOptionsGetReal(NULL, NULL, "-jump_magnitude", &opts_mag, NULL);
1985:           DMDAGetInfo(da_Stokes, NULL, NULL, NULL, NULL, &px, &py, &pz, NULL, NULL, NULL, NULL, NULL, NULL);
1986:           rr = PetscGlobalRank % (px * py);
1987:           if (px % 2) jump = (PetscBool)(rr % 2);
1988:           else jump = (PetscBool)((rr / px) % 2 ? rr % 2 : !(rr % 2));
1989:           rr = PetscGlobalRank / (px * py);
1990:           if (rr % 2) jump = (PetscBool)!jump;
1991:           CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell);
1992:           for (p = 0; p < GAUSS_POINTS; p++) {
1993:             PetscReal xp = PetscRealPart(cell->gp_coords[NSD * p]);
1994:             PetscReal yp = PetscRealPart(cell->gp_coords[NSD * p + 1]);

1996:             cell->eta[p] = jump ? PetscPowReal(10.0, opts_mag) : opts_eta0;
1997:             cell->fx[p]  = 0.0;
1998:             cell->fy[p]  = -PetscSinReal(2.2 * PETSC_PI * yp) * PetscCosReal(1.0 * PETSC_PI * xp);
1999:             cell->fz[p]  = 0.0;
2000:             cell->hc[p]  = 0.0;
2001:           }
2002:         }
2003:       }
2004:     }
2005:     break;
2006:   default:
2007:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "No default model is supported. Choose either -model {0,1,2,3}");
2008:   }

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

2012:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
2013:   DMCreateMatrix(da_Stokes, &A);
2014:   DMCreateMatrix(da_Stokes, &B);
2015:   DMCreateGlobalVector(da_Stokes, &X);
2016:   DMCreateGlobalVector(da_Stokes, &f);

2018:   /* assemble A11 */
2019:   MatZeroEntries(A);
2020:   MatZeroEntries(B);
2021:   VecZeroEntries(f);

2023:   AssembleA_Stokes(A, da_Stokes, cell_properties);
2024:   MatViewFromOptions(A, NULL, "-amat_view");
2025:   AssembleA_PCStokes(B, da_Stokes, cell_properties);
2026:   MatViewFromOptions(B, NULL, "-bmat_view");
2027:   /* build force vector */
2028:   AssembleF_Stokes(f, da_Stokes, cell_properties);

2030:   /* SOLVE */
2031:   KSPCreate(PETSC_COMM_WORLD, &ksp_S);
2032:   KSPSetOptionsPrefix(ksp_S, "stokes_");
2033:   KSPSetOperators(ksp_S, A, B);
2034:   KSPSetFromOptions(ksp_S);

2036:   {
2037:     PC             pc;
2038:     const PetscInt ufields[] = {0, 1, 2}, pfields[] = {3};
2039:     KSPGetPC(ksp_S, &pc);
2040:     PCFieldSplitSetBlockSize(pc, 4);
2041:     PCFieldSplitSetFields(pc, "u", 3, ufields, ufields);
2042:     PCFieldSplitSetFields(pc, "p", 1, pfields, pfields);
2043:   }

2045:   {
2046:     PC        pc;
2047:     PetscBool same = PETSC_FALSE;
2048:     KSPGetPC(ksp_S, &pc);
2049:     PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
2050:     if (same) PCMGSetupViaCoarsen(pc, da_Stokes);
2051:   }

2053:   {
2054:     PC        pc;
2055:     PetscBool same = PETSC_FALSE;
2056:     KSPGetPC(ksp_S, &pc);
2057:     PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &same);
2058:     if (same) KSPSetOperators(ksp_S, A, A);
2059:   }

2061:   {
2062:     PetscBool stokes_monitor = PETSC_FALSE;
2063:     PetscOptionsGetBool(NULL, NULL, "-stokes_ksp_monitor_blocks", &stokes_monitor, 0);
2064:     if (stokes_monitor) KSPMonitorSet(ksp_S, KSPMonitorStokesBlocks, NULL, NULL);
2065:   }

2067:   if (resolve) {
2068:     /* Test changing matrix data structure and resolve */
2069:     VecDuplicate(X, &X1);
2070:   }

2072:   KSPSolve(ksp_S, f, X);
2073:   if (resolve) {
2074:     Mat C;
2075:     MatDuplicate(A, MAT_COPY_VALUES, &C);
2076:     KSPSetOperators(ksp_S, C, C);
2077:     KSPSolve(ksp_S, f, X1);
2078:     MatDestroy(&C);
2079:     VecDestroy(&X1);
2080:   }

2082:   PetscOptionsGetBool(NULL, NULL, "-write_pvts", &write_output, NULL);
2083:   if (write_output) DAView3DPVTS(da_Stokes, X, "up");
2084:   {
2085:     PetscBool flg = PETSC_FALSE;
2086:     char      filename[PETSC_MAX_PATH_LEN];
2087:     PetscOptionsGetString(NULL, NULL, "-write_binary", filename, sizeof(filename), &flg);
2088:     if (flg) {
2089:       PetscViewer viewer;
2090:       /* PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer); */
2091:       PetscViewerVTKOpen(PETSC_COMM_WORLD, "ex42.vts", FILE_MODE_WRITE, &viewer);
2092:       VecView(X, viewer);
2093:       PetscViewerDestroy(&viewer);
2094:     }
2095:   }
2096:   KSPGetIterationNumber(ksp_S, &its);

2098:   /* verify */
2099:   if (model_definition == 1) {
2100:     DM  da_Stokes_analytic;
2101:     Vec X_analytic;

2103:     DMDACreateManufacturedSolution(mx, my, mz, &da_Stokes_analytic, &X_analytic);
2104:     if (write_output) DAView3DPVTS(da_Stokes_analytic, X_analytic, "ms");
2105:     DMDAIntegrateErrors3D(da_Stokes_analytic, X, X_analytic);
2106:     if (write_output) DAView3DPVTS(da_Stokes, X, "up2");
2107:     DMDestroy(&da_Stokes_analytic);
2108:     VecDestroy(&X_analytic);
2109:   }

2111:   KSPDestroy(&ksp_S);
2112:   VecDestroy(&X);
2113:   VecDestroy(&f);
2114:   MatDestroy(&A);
2115:   MatDestroy(&B);

2117:   CellPropertiesDestroy(&cell_properties);
2118:   DMDestroy(&da_Stokes);
2119:   return 0;
2120: }

2122: int main(int argc, char **args)
2123: {
2124:   PetscInt mx, my, mz;

2127:   PetscInitialize(&argc, &args, (char *)0, help);
2128:   mx = my = mz = 10;
2129:   PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL);
2130:   my = mx;
2131:   mz = mx;
2132:   PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL);
2133:   PetscOptionsGetInt(NULL, NULL, "-mz", &mz, NULL);
2134:   solve_stokes_3d_coupled(mx, my, mz);
2135:   PetscFinalize();
2136:   return 0;
2137: }

2139: /*TEST

2141:    test:
2142:       args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type lu

2144:    test:
2145:       suffix: 2
2146:       nsize: 3
2147:       args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type redundant -stokes_redundant_pc_type lu

2149:    test:
2150:       suffix: bddc_stokes
2151:       nsize: 6
2152:       args: -mx 5 -my 4 -mz 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_dirichlet_pc_type svd -stokes_pc_bddc_neumann_pc_type svd -stokes_pc_bddc_coarse_redundant_pc_type svd

2154:    test:
2155:       suffix: bddc_stokes_deluxe
2156:       nsize: 6
2157:       args: -mx 5 -my 4 -mz 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_dirichlet_pc_type svd -stokes_pc_bddc_neumann_pc_type svd -stokes_pc_bddc_coarse_redundant_pc_type svd -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc

2159:    test:
2160:       requires: !single
2161:       suffix: bddc_stokes_subdomainjump_deluxe
2162:       nsize: 9
2163:       args: -model 4 -jump_magnitude 4 -mx 6 -my 6 -mz 2 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc -stokes_pc_bddc_schur_layers 1

2165:    test:
2166:       requires: !single
2167:       suffix: 3
2168:       args: -stokes_ksp_converged_reason -stokes_pc_type fieldsplit -resolve

2170:    test:
2171:       suffix: tut_1
2172:       nsize: 4
2173:       requires: !single
2174:       args: -stokes_ksp_monitor

2176:    test:
2177:       suffix: tut_2
2178:       nsize: 4
2179:       requires: !single
2180:       args: -stokes_ksp_monitor -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_type schur

2182:    test:
2183:       suffix: tut_3
2184:       nsize: 4
2185:       requires: !single
2186:       args: -mx 20 -stokes_ksp_monitor -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_type schur

2188: TEST*/