Actual source code: ex42.c

petsc-3.10.2 2018-10-09
Report Typos and Errors
  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: {
 53:   CellProperties cells;
 54:   PetscInt       mx,my,mz,sex,sey,sez;

 57:   PetscNew(&cells);

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

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

 71:   *C = cells;
 72:   return(0);
 73: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

209:   JJ[0][0] = J00;      JJ[0][1] = J01;      JJ[0][2] = J02;
210:   JJ[1][0] = J10;      JJ[1][1] = J11;      JJ[1][2] = J12;
211:   JJ[2][0] = J20;      JJ[2][1] = J21;      JJ[2][2] = J22;

213:   matrix_inverse_3x3(JJ,iJ);

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

217:   for (n=0; n<NODES_PER_EL; n++) {
218:     GNx[0][n] = GNi[0][n]*iJ[0][0] + GNi[1][n]*iJ[0][1] + GNi[2][n]*iJ[0][2];
219:     GNx[1][n] = GNi[0][n]*iJ[1][0] + GNi[1][n]*iJ[1][1] + GNi[2][n]*iJ[1][2];
220:     GNx[2][n] = GNi[0][n]*iJ[2][0] + GNi[1][n]*iJ[2][1] + GNi[2][n]*iJ[2][2];
221:   }
222: }

224: static void ConstructGaussQuadrature3D(PetscInt *ngp,PetscScalar gp_xi[][NSD],PetscScalar gp_weight[])
225: {
226:   *ngp        = 8;
227:   gp_xi[0][0] = -0.57735026919; gp_xi[0][1] = -0.57735026919; gp_xi[0][2] = -0.57735026919;
228:   gp_xi[1][0] = -0.57735026919; gp_xi[1][1] =  0.57735026919; gp_xi[1][2] = -0.57735026919;
229:   gp_xi[2][0] =  0.57735026919; gp_xi[2][1] =  0.57735026919; gp_xi[2][2] = -0.57735026919;
230:   gp_xi[3][0] =  0.57735026919; gp_xi[3][1] = -0.57735026919; gp_xi[3][2] = -0.57735026919;

232:   gp_xi[4][0] = -0.57735026919; gp_xi[4][1] = -0.57735026919; gp_xi[4][2] =  0.57735026919;
233:   gp_xi[5][0] = -0.57735026919; gp_xi[5][1] =  0.57735026919; gp_xi[5][2] =  0.57735026919;
234:   gp_xi[6][0] =  0.57735026919; gp_xi[6][1] =  0.57735026919; gp_xi[6][2] =  0.57735026919;
235:   gp_xi[7][0] =  0.57735026919; gp_xi[7][1] = -0.57735026919; gp_xi[7][2] =  0.57735026919;

237:   gp_weight[0] = 1.0;
238:   gp_weight[1] = 1.0;
239:   gp_weight[2] = 1.0;
240:   gp_weight[3] = 1.0;

242:   gp_weight[4] = 1.0;
243:   gp_weight[5] = 1.0;
244:   gp_weight[6] = 1.0;
245:   gp_weight[7] = 1.0;
246: }

248: /*
249:  i,j are the element indices
250:  The unknown is a vector quantity.
251:  The s[].c is used to indicate the degree of freedom.
252:  */
253: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j,PetscInt k)
254: {
255:   PetscInt n;

258:   /* velocity */
259:   n = 0;
260:   /* node 0 */
261:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++; /* Vx0 */
262:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++; /* Vy0 */
263:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++; /* Vz0 */

265:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
266:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
267:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;

269:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
270:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
271:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;

273:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++;
274:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++;
275:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++;

277:   /* */
278:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++; /* Vx4 */
279:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++; /* Vy4 */
280:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++; /* Vz4 */

282:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
283:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
284:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;

286:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
287:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
288:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;

290:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++;
291:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++;
292:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++;

294:   /* pressure */
295:   n = 0;

297:   s_p[n].i = i;   s_p[n].j = j;   s_p[n].k = k; s_p[n].c = 3; n++; /* P0 */
298:   s_p[n].i = i;   s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
299:   s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
300:   s_p[n].i = i+1; s_p[n].j = j;   s_p[n].k = k; s_p[n].c = 3; n++;

302:   s_p[n].i = i;   s_p[n].j = j;   s_p[n].k = k+1; s_p[n].c = 3; n++; /* P0 */
303:   s_p[n].i = i;   s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
304:   s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
305:   s_p[n].i = i+1; s_p[n].j = j;   s_p[n].k = k+1; s_p[n].c = 3; n++;
306:   return(0);
307: }

309: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords,PetscInt i,PetscInt j,PetscInt k,PetscScalar el_coord[])
310: {
312:   /* get coords for the element */
313:   el_coord[0] = coords[k][j][i].x;
314:   el_coord[1] = coords[k][j][i].y;
315:   el_coord[2] = coords[k][j][i].z;

317:   el_coord[3] = coords[k][j+1][i].x;
318:   el_coord[4] = coords[k][j+1][i].y;
319:   el_coord[5] = coords[k][j+1][i].z;

321:   el_coord[6] = coords[k][j+1][i+1].x;
322:   el_coord[7] = coords[k][j+1][i+1].y;
323:   el_coord[8] = coords[k][j+1][i+1].z;

325:   el_coord[9]  = coords[k][j][i+1].x;
326:   el_coord[10] = coords[k][j][i+1].y;
327:   el_coord[11] = coords[k][j][i+1].z;

329:   el_coord[12] = coords[k+1][j][i].x;
330:   el_coord[13] = coords[k+1][j][i].y;
331:   el_coord[14] = coords[k+1][j][i].z;

333:   el_coord[15] = coords[k+1][j+1][i].x;
334:   el_coord[16] = coords[k+1][j+1][i].y;
335:   el_coord[17] = coords[k+1][j+1][i].z;

337:   el_coord[18] = coords[k+1][j+1][i+1].x;
338:   el_coord[19] = coords[k+1][j+1][i+1].y;
339:   el_coord[20] = coords[k+1][j+1][i+1].z;

341:   el_coord[21] = coords[k+1][j][i+1].x;
342:   el_coord[22] = coords[k+1][j][i+1].y;
343:   el_coord[23] = coords[k+1][j][i+1].z;
344:   return(0);
345: }

347: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field,PetscInt i,PetscInt j,PetscInt k,StokesDOF nodal_fields[])
348: {
350:   /* get the nodal fields for u */
351:   nodal_fields[0].u_dof = field[k][j][i].u_dof;
352:   nodal_fields[0].v_dof = field[k][j][i].v_dof;
353:   nodal_fields[0].w_dof = field[k][j][i].w_dof;

355:   nodal_fields[1].u_dof = field[k][j+1][i].u_dof;
356:   nodal_fields[1].v_dof = field[k][j+1][i].v_dof;
357:   nodal_fields[1].w_dof = field[k][j+1][i].w_dof;

359:   nodal_fields[2].u_dof = field[k][j+1][i+1].u_dof;
360:   nodal_fields[2].v_dof = field[k][j+1][i+1].v_dof;
361:   nodal_fields[2].w_dof = field[k][j+1][i+1].w_dof;

363:   nodal_fields[3].u_dof = field[k][j][i+1].u_dof;
364:   nodal_fields[3].v_dof = field[k][j][i+1].v_dof;
365:   nodal_fields[3].w_dof = field[k][j][i+1].w_dof;

367:   nodal_fields[4].u_dof = field[k+1][j][i].u_dof;
368:   nodal_fields[4].v_dof = field[k+1][j][i].v_dof;
369:   nodal_fields[4].w_dof = field[k+1][j][i].w_dof;

371:   nodal_fields[5].u_dof = field[k+1][j+1][i].u_dof;
372:   nodal_fields[5].v_dof = field[k+1][j+1][i].v_dof;
373:   nodal_fields[5].w_dof = field[k+1][j+1][i].w_dof;

375:   nodal_fields[6].u_dof = field[k+1][j+1][i+1].u_dof;
376:   nodal_fields[6].v_dof = field[k+1][j+1][i+1].v_dof;
377:   nodal_fields[6].w_dof = field[k+1][j+1][i+1].w_dof;

379:   nodal_fields[7].u_dof = field[k+1][j][i+1].u_dof;
380:   nodal_fields[7].v_dof = field[k+1][j][i+1].v_dof;
381:   nodal_fields[7].w_dof = field[k+1][j][i+1].w_dof;

383:   /* get the nodal fields for p */
384:   nodal_fields[0].p_dof = field[k][j][i].p_dof;
385:   nodal_fields[1].p_dof = field[k][j+1][i].p_dof;
386:   nodal_fields[2].p_dof = field[k][j+1][i+1].p_dof;
387:   nodal_fields[3].p_dof = field[k][j][i+1].p_dof;

389:   nodal_fields[4].p_dof = field[k+1][j][i].p_dof;
390:   nodal_fields[5].p_dof = field[k+1][j+1][i].p_dof;
391:   nodal_fields[6].p_dof = field[k+1][j+1][i+1].p_dof;
392:   nodal_fields[7].p_dof = field[k+1][j][i+1].p_dof;
393:   return(0);
394: }

396: 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)
397: {
398:   PetscInt              ij;
399:   PETSC_UNUSED PetscInt r,c,nr,nc;

401:   nr = w_NPE*w_dof;
402:   nc = u_NPE*u_dof;

404:   r = w_dof*wi+wd;
405:   c = u_dof*ui+ud;

407:   ij = r*nc+c;

409:   return ij;
410: }

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

417:   for (n = 0; n<NODES_PER_EL; n++) {
418:     II = u_eqn[NSD*n].i;
419:     J = u_eqn[NSD*n].j;
420:     K = u_eqn[NSD*n].k;

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

424:     II = u_eqn[NSD*n+1].i;
425:     J = u_eqn[NSD*n+1].j;
426:     K = u_eqn[NSD*n+1].k;

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

430:     II = u_eqn[NSD*n+2].i;
431:     J = u_eqn[NSD*n+2].j;
432:     K = u_eqn[NSD*n+2].k;
433:     fields_F[K][J][II].w_dof = fields_F[K][J][II].w_dof+Fe_u[NSD*n+2];

435:     II = p_eqn[n].i;
436:     J = p_eqn[n].j;
437:     K = p_eqn[n].k;

439:     fields_F[K][J][II].p_dof = fields_F[K][J][II].p_dof+Fe_p[n];

441:   }
442:   return(0);
443: }

445: static void FormStressOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
446: {
447:   PetscInt       ngp;
448:   PetscScalar    gp_xi[GAUSS_POINTS][NSD];
449:   PetscScalar    gp_weight[GAUSS_POINTS];
450:   PetscInt       p,i,j,k;
451:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
452:   PetscScalar    J_p,tildeD[6];
453:   PetscScalar    B[6][U_DOFS*NODES_PER_EL];
454:   const PetscInt nvdof = U_DOFS*NODES_PER_EL;

456:   /* define quadrature rule */
457:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

459:   /* evaluate integral */
460:   for (p = 0; p < ngp; p++) {
461:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
462:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);

464:     for (i = 0; i < NODES_PER_EL; i++) {
465:       PetscScalar d_dx_i = GNx_p[0][i];
466:       PetscScalar d_dy_i = GNx_p[1][i];
467:       PetscScalar d_dz_i = GNx_p[2][i];

469:       B[0][3*i] = d_dx_i; B[0][3*i+1] = 0.0;     B[0][3*i+2] = 0.0;
470:       B[1][3*i] = 0.0;    B[1][3*i+1] = d_dy_i;  B[1][3*i+2] = 0.0;
471:       B[2][3*i] = 0.0;    B[2][3*i+1] = 0.0;     B[2][3*i+2] = d_dz_i;

473:       B[3][3*i] = d_dy_i; B[3][3*i+1] = d_dx_i;  B[3][3*i+2] = 0.0;   /* e_xy */
474:       B[4][3*i] = d_dz_i; B[4][3*i+1] = 0.0;     B[4][3*i+2] = d_dx_i; /* e_xz */
475:       B[5][3*i] = 0.0;    B[5][3*i+1] = d_dz_i;  B[5][3*i+2] = d_dy_i; /* e_yz */
476:     }


479:     tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
480:     tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
481:     tildeD[2] = 2.0*gp_weight[p]*J_p*eta[p];

483:     tildeD[3] =     gp_weight[p]*J_p*eta[p];
484:     tildeD[4] =     gp_weight[p]*J_p*eta[p];
485:     tildeD[5] =     gp_weight[p]*J_p*eta[p];

487:     /* form Bt tildeD B */
488:     /*
489:      Ke_ij = Bt_ik . D_kl . B_lj
490:      = B_ki . D_kl . B_lj
491:      = B_ki . D_kk . B_kj
492:      */
493:     for (i = 0; i < nvdof; i++) {
494:       for (j = i; j < nvdof; j++) {
495:         for (k = 0; k < 6; k++) {
496:           Ke[i*nvdof+j] += B[k][i]*tildeD[k]*B[k][j];
497:         }
498:       }
499:     }

501:   }
502:   /* fill lower triangular part */
503: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
504:   for (i = 0; i < nvdof; i++) {
505:     for (j = i; j < nvdof; j++) {
506:       Ke[j*nvdof+i] = Ke[i*nvdof+j];
507:     }
508:   }
509: #endif
510: }

512: static void FormGradientOperatorQ13D(PetscScalar Ke[],PetscScalar coords[])
513: {
514:   PetscInt    ngp;
515:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
516:   PetscScalar gp_weight[GAUSS_POINTS];
517:   PetscInt    p,i,j,di;
518:   PetscScalar Ni_p[NODES_PER_EL];
519:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
520:   PetscScalar J_p,fac;

522:   /* define quadrature rule */
523:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

525:   /* evaluate integral */
526:   for (p = 0; p < ngp; p++) {
527:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
528:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
529:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
530:     fac = gp_weight[p]*J_p;

532:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
533:       for (di = 0; di < NSD; di++) { /* u dofs */
534:         for (j = 0; j < NODES_PER_EL; j++) {  /* p nodes, p dofs = 1 (ie no loop) */
535:           PetscInt IJ;
536:           IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,3,j,0,NODES_PER_EL,1);

538:           Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
539:         }
540:       }
541:     }
542:   }
543: }

545: static void FormDivergenceOperatorQ13D(PetscScalar De[],PetscScalar coords[])
546: {
547:   PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
548:   PetscInt    i,j;
549:   PetscInt    nr_g,nc_g;

551:   PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
552:   FormGradientOperatorQ13D(Ge,coords);

554:   nr_g = U_DOFS*NODES_PER_EL;
555:   nc_g = P_DOFS*NODES_PER_EL;

557:   for (i = 0; i < nr_g; i++) {
558:     for (j = 0; j < nc_g; j++) {
559:       De[nr_g*j+i] = Ge[nc_g*i+j];
560:     }
561:   }
562: }

564: static void FormStabilisationOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
565: {
566:   PetscInt    ngp;
567:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
568:   PetscScalar gp_weight[GAUSS_POINTS];
569:   PetscInt    p,i,j;
570:   PetscScalar Ni_p[NODES_PER_EL];
571:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
572:   PetscScalar J_p,fac,eta_avg;

574:   /* define quadrature rule */
575:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

577:   /* evaluate integral */
578:   for (p = 0; p < ngp; p++) {
579:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
580:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
581:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
582:     fac = gp_weight[p]*J_p;
583:     /*
584:      for (i = 0; i < NODES_PER_EL; i++) {
585:      for (j = i; j < NODES_PER_EL; j++) {
586:      Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
587:      }
588:      }
589:      */

591:     for (i = 0; i < NODES_PER_EL; i++) {
592:       for (j = 0; j < NODES_PER_EL; j++) {
593:         Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
594:       }
595:     }
596:   }

598:   /* scale */
599:   eta_avg = 0.0;
600:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
601:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
602:   fac     = 1.0/eta_avg;

604:   /*
605:    for (i = 0; i < NODES_PER_EL; i++) {
606:    for (j = i; j < NODES_PER_EL; j++) {
607:    Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
608:    #if defined(ASSEMBLE_LOWER_TRIANGULAR)
609:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
610:    #endif
611:    }
612:    }
613:    */

615:   for (i = 0; i < NODES_PER_EL; i++) {
616:     for (j = 0; j < NODES_PER_EL; j++) {
617:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
618:     }
619:   }
620: }

622: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
623: {
624:   PetscInt    ngp;
625:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
626:   PetscScalar gp_weight[GAUSS_POINTS];
627:   PetscInt    p,i,j;
628:   PetscScalar Ni_p[NODES_PER_EL];
629:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
630:   PetscScalar J_p,fac,eta_avg;

632:   /* define quadrature rule */
633:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

635:   /* evaluate integral */
636:   for (p = 0; p < ngp; p++) {
637:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
638:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
639:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
640:     fac = gp_weight[p]*J_p;

642:     /*
643:      for (i = 0; i < NODES_PER_EL; i++) {
644:      for (j = i; j < NODES_PER_EL; j++) {
645:      Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
646:      }
647:      }
648:      */

650:     for (i = 0; i < NODES_PER_EL; i++) {
651:       for (j = 0; j < NODES_PER_EL; j++) {
652:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
653:       }
654:     }
655:   }

657:   /* scale */
658:   eta_avg = 0.0;
659:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
660:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
661:   fac     = 1.0/eta_avg;
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*Ke[NODES_PER_EL*i+j];
666:    #if defined(ASSEMBLE_LOWER_TRIANGULAR)
667:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
668:    #endif
669:    }
670:    }
671:    */

673:   for (i = 0; i < NODES_PER_EL; i++) {
674:     for (j = 0; j < NODES_PER_EL; j++) {
675:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
676:     }
677:   }
678: }

680: static void FormMomentumRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[],PetscScalar fz[])
681: {
682:   PetscInt    ngp;
683:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
684:   PetscScalar gp_weight[GAUSS_POINTS];
685:   PetscInt    p,i;
686:   PetscScalar Ni_p[NODES_PER_EL];
687:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
688:   PetscScalar J_p,fac;

690:   /* define quadrature rule */
691:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

693:   /* evaluate integral */
694:   for (p = 0; p < ngp; p++) {
695:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
696:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
697:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
698:     fac = gp_weight[p]*J_p;

700:     for (i = 0; i < NODES_PER_EL; i++) {
701:       Fe[NSD*i]   -= fac*Ni_p[i]*fx[p];
702:       Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
703:       Fe[NSD*i+2] -= fac*Ni_p[i]*fz[p];
704:     }
705:   }
706: }

708: static void FormContinuityRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar hc[])
709: {
710:   PetscInt    ngp;
711:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
712:   PetscScalar gp_weight[GAUSS_POINTS];
713:   PetscInt    p,i;
714:   PetscScalar Ni_p[NODES_PER_EL];
715:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
716:   PetscScalar J_p,fac;

718:   /* define quadrature rule */
719:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

721:   /* evaluate integral */
722:   for (p = 0; p < ngp; p++) {
723:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
724:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
725:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
726:     fac = gp_weight[p]*J_p;

728:     for (i = 0; i < NODES_PER_EL; i++) Fe[i] -= fac*Ni_p[i]*hc[p];
729:   }
730: }

732: #define _ZERO_ROWCOL_i(A,i) {                   \
733:     PetscInt    KK;                             \
734:     PetscScalar tmp = A[24*(i)+(i)];            \
735:     for (KK=0;KK<24;KK++) A[24*(i)+KK]=0.0;     \
736:     for (KK=0;KK<24;KK++) A[24*KK+(i)]=0.0;     \
737:     A[24*(i)+(i)] = tmp;}                       \

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

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

747: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,CellProperties cell_properties)
748: {
749:   DM                     cda;
750:   Vec                    coords;
751:   DMDACoor3d             ***_coords;
752:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS];
753:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS];
754:   PetscInt               sex,sey,sez,mx,my,mz;
755:   PetscInt               ei,ej,ek;
756:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
757:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
758:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
759:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
760:   PetscScalar            el_coords[NODES_PER_EL*NSD];
761:   GaussPointCoefficients *props;
762:   PetscScalar            *prop_eta;
763:   PetscInt               n,M,N,P;
764:   PetscErrorCode         ierr;

767:   DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
768:   /* setup for coords */
769:   DMGetCoordinateDM(stokes_da,&cda);
770:   DMGetCoordinatesLocal(stokes_da,&coords);
771:   DMDAVecGetArray(cda,coords,&_coords);
772:   DMDAGetElementsCorners(stokes_da,&sex,&sey,&sez);
773:   DMDAGetElementsSizes(stokes_da,&mx,&my,&mz);
774:   for (ek = sez; ek < sez+mz; ek++) {
775:     for (ej = sey; ej < sey+my; ej++) {
776:       for (ei = sex; ei < sex+mx; ei++) {
777:         /* get coords for the element */
778:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
779:         /* get cell properties */
780:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
781:         /* get coefficients for the element */
782:         prop_eta = props->eta;

784:         /* initialise element stiffness matrix */
785:         PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
786:         PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
787:         PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
788:         PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);

790:         /* form element stiffness matrix */
791:         FormStressOperatorQ13D(Ae,el_coords,prop_eta);
792:         FormGradientOperatorQ13D(Ge,el_coords);
793:         /*#if defined(ASSEMBLE_LOWER_TRIANGULAR)*/
794:         FormDivergenceOperatorQ13D(De,el_coords);
795:         /*#endif*/
796:         FormStabilisationOperatorQ13D(Ce,el_coords,prop_eta);

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

801:         for (n=0; n<NODES_PER_EL; n++) {
802:           if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
803:             _ZERO_ROWCOL_i(Ae,3*n);
804:             _ZERO_ROW_i   (Ge,3*n);
805:             _ZERO_COL_i   (De,3*n);
806:           }

808:           if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
809:             _ZERO_ROWCOL_i(Ae,3*n+1);
810:             _ZERO_ROW_i   (Ge,3*n+1);
811:             _ZERO_COL_i   (De,3*n+1);
812:           }

814:           if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
815:             _ZERO_ROWCOL_i(Ae,3*n+2);
816:             _ZERO_ROW_i   (Ge,3*n+2);
817:             _ZERO_COL_i   (De,3*n+2);
818:           }
819:         }
820:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
821:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
822:         MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
823:         MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
824:       }
825:     }
826:   }
827:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
828:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

832:   return(0);
833: }

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

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

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

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

879:         /* form element stiffness matrix */
880:         FormStressOperatorQ13D(Ae,el_coords,prop_eta);
881:         FormGradientOperatorQ13D(Ge,el_coords);
882:         /* FormDivergenceOperatorQ13D(De,el_coords); */
883:         FormScaledMassMatrixOperatorQ13D(Ce,el_coords,prop_eta);

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

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

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

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

917:   DMDAVecRestoreArray(cda,coords,&_coords);
918:   return(0);
919: }

921: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,CellProperties cell_properties)
922: {
923:   DM                     cda;
924:   Vec                    coords;
925:   DMDACoor3d             ***_coords;
926:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS];
927:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS];
928:   PetscInt               sex,sey,sez,mx,my,mz;
929:   PetscInt               ei,ej,ek;
930:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
931:   PetscScalar            He[NODES_PER_EL*P_DOFS];
932:   PetscScalar            el_coords[NODES_PER_EL*NSD];
933:   GaussPointCoefficients *props;
934:   PetscScalar            *prop_fx,*prop_fy,*prop_fz,*prop_hc;
935:   Vec                    local_F;
936:   StokesDOF              ***ff;
937:   PetscInt               n,M,N,P;
938:   PetscErrorCode         ierr;

941:   DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
942:   /* setup for coords */
943:   DMGetCoordinateDM(stokes_da,&cda);
944:   DMGetCoordinatesLocal(stokes_da,&coords);
945:   DMDAVecGetArray(cda,coords,&_coords);

947:   /* get acces to the vector */
948:   DMGetLocalVector(stokes_da,&local_F);
949:   VecZeroEntries(local_F);
950:   DMDAVecGetArray(stokes_da,local_F,&ff);
951:   DMDAGetElementsCorners(stokes_da,&sex,&sey,&sez);
952:   DMDAGetElementsSizes(stokes_da,&mx,&my,&mz);
953:   for (ek = sez; ek < sez+mz; ek++) {
954:     for (ej = sey; ej < sey+my; ej++) {
955:       for (ei = sex; ei < sex+mx; ei++) {
956:         /* get coords for the element */
957:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
958:         /* get cell properties */
959:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
960:         /* get coefficients for the element */
961:         prop_fx = props->fx;
962:         prop_fy = props->fy;
963:         prop_fz = props->fz;
964:         prop_hc = props->hc;

966:         /* initialise element stiffness matrix */
967:         PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
968:         PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);

970:         /* form element stiffness matrix */
971:         FormMomentumRhsQ13D(Fe,el_coords,prop_fx,prop_fy,prop_fz);
972:         FormContinuityRhsQ13D(He,el_coords,prop_hc);

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

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

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

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

985:         DMDASetValuesLocalStencil3D_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
986:       }
987:     }
988:   }
989:   DMDAVecRestoreArray(stokes_da,local_F,&ff);
990:   DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
991:   DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
992:   DMRestoreLocalVector(stokes_da,&local_F);

994:   DMDAVecRestoreArray(cda,coords,&_coords);
995:   return(0);
996: }

998: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta,PetscReal *MX,PetscReal *MY,PetscReal *MZ)
999: {
1000:   *theta = 0.0;
1001:   *MX    = 2.0 * PETSC_PI;
1002:   *MY    = 2.0 * PETSC_PI;
1003:   *MZ    = 2.0 * PETSC_PI;
1004: }
1005: static void evaluate_MS_FrankKamentski(PetscReal pos[],PetscReal v[],PetscReal *p,PetscReal *eta,PetscReal Fm[],PetscReal *Fc)
1006: {
1007:   PetscReal x,y,z;
1008:   PetscReal theta,MX,MY,MZ;

1010:   evaluate_MS_FrankKamentski_constants(&theta,&MX,&MY,&MZ);
1011:   x = pos[0];
1012:   y = pos[1];
1013:   z = pos[2];
1014:   if (v) {
1015:     /*
1016:      v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1017:      v[1] = z*cos(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1018:      v[2] = -(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z);
1019:      */
1020:     /*
1021:      v[0] = PetscSinReal(PETSC_PI*x);
1022:      v[1] = PetscSinReal(PETSC_PI*y);
1023:      v[2] = PetscSinReal(PETSC_PI*z);
1024:      */
1025:     v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1026:     v[1] = z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1027:     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;
1028:   }
1029:   if (p) *p = PetscPowRealInt(x,2) + PetscPowRealInt(y,2) + PetscPowRealInt(z,2);
1030:   if (eta) {
1031:     /**eta = PetscExpReal(-theta*(1.0 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));*/
1032:     *eta = 1.0;
1033:   }
1034:   if (Fm) {
1035:     /*
1036:      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))) ;
1037:      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)));
1038:      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)))  ;
1039:      */
1040:     /*
1041:      Fm[0]=-2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*x);
1042:      Fm[1]=-2*y - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*y);
1043:      Fm[2]=-2*z - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*z);
1044:      */
1045:     /*
1046:      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) ;
1047:      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);
1048:      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) ;
1049:      */
1050:     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);
1051:     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);
1052:     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);
1053:   }
1054:   if (Fc) {
1055:     /**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) ;*/
1056:     /**Fc = PETSC_PI*PetscCosReal(PETSC_PI*x) + PETSC_PI*PetscCosReal(PETSC_PI*y) + PETSC_PI*PetscCosReal(PETSC_PI*z);*/
1057:     /**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);*/
1058:     *Fc = 0.0;
1059:   }
1060: }

1062: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx,PetscInt my,PetscInt mz,DM *_da,Vec *_X)
1063: {
1064:   DM             da,cda;
1065:   Vec            X;
1066:   StokesDOF      ***_stokes;
1067:   Vec            coords;
1068:   DMDACoor3d     ***_coords;
1069:   PetscInt       si,sj,sk,ei,ej,ek,i,j,k;

1073:   DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1074:                       mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,4,1,NULL,NULL,NULL,&da);
1075:   DMSetFromOptions(da);
1076:   DMSetUp(da);
1077:   DMDASetFieldName(da,0,"anlytic_Vx");
1078:   DMDASetFieldName(da,1,"anlytic_Vy");
1079:   DMDASetFieldName(da,2,"anlytic_Vz");
1080:   DMDASetFieldName(da,3,"analytic_P");

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

1084:   DMGetCoordinatesLocal(da,&coords);
1085:   DMGetCoordinateDM(da,&cda);
1086:   DMDAVecGetArray(cda,coords,&_coords);

1088:   DMCreateGlobalVector(da,&X);
1089:   DMDAVecGetArray(da,X,&_stokes);

1091:   DMDAGetCorners(da,&si,&sj,&sk,&ei,&ej,&ek);
1092:   for (k = sk; k < sk+ek; k++) {
1093:     for (j = sj; j < sj+ej; j++) {
1094:       for (i = si; i < si+ei; i++) {
1095:         PetscReal pos[NSD],pressure,vel[NSD];

1097:         pos[0] = PetscRealPart(_coords[k][j][i].x);
1098:         pos[1] = PetscRealPart(_coords[k][j][i].y);
1099:         pos[2] = PetscRealPart(_coords[k][j][i].z);

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

1103:         _stokes[k][j][i].u_dof = vel[0];
1104:         _stokes[k][j][i].v_dof = vel[1];
1105:         _stokes[k][j][i].w_dof = vel[2];
1106:         _stokes[k][j][i].p_dof = pressure;
1107:       }
1108:     }
1109:   }
1110:   DMDAVecRestoreArray(da,X,&_stokes);
1111:   DMDAVecRestoreArray(cda,coords,&_coords);

1113:   *_da = da;
1114:   *_X  = X;
1115:   return(0);
1116: }

1118: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da,Vec X,Vec X_analytic)
1119: {
1120:   DM          cda;
1121:   Vec         coords,X_analytic_local,X_local;
1122:   DMDACoor3d  ***_coords;
1123:   PetscInt    sex,sey,sez,mx,my,mz;
1124:   PetscInt    ei,ej,ek;
1125:   PetscScalar el_coords[NODES_PER_EL*NSD];
1126:   StokesDOF   ***stokes_analytic,***stokes;
1127:   StokesDOF   stokes_analytic_e[NODES_PER_EL],stokes_e[NODES_PER_EL];

1129:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1130:   PetscScalar    Ni_p[NODES_PER_EL];
1131:   PetscInt       ngp;
1132:   PetscScalar    gp_xi[GAUSS_POINTS][NSD];
1133:   PetscScalar    gp_weight[GAUSS_POINTS];
1134:   PetscInt       p,i;
1135:   PetscScalar    J_p,fac;
1136:   PetscScalar    h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1137:   PetscScalar    tint_p_ms,tint_p,int_p_ms,int_p;
1138:   PetscInt       M;
1139:   PetscReal      xymin[NSD],xymax[NSD];

1143:   /* define quadrature rule */
1144:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

1146:   /* setup for coords */
1147:   DMGetCoordinateDM(stokes_da,&cda);
1148:   DMGetCoordinatesLocal(stokes_da,&coords);
1149:   DMDAVecGetArray(cda,coords,&_coords);

1151:   /* setup for analytic */
1152:   DMCreateLocalVector(stokes_da,&X_analytic_local);
1153:   DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1154:   DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1155:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1157:   /* setup for solution */
1158:   DMCreateLocalVector(stokes_da,&X_local);
1159:   DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1160:   DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1161:   DMDAVecGetArray(stokes_da,X_local,&stokes);

1163:   DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1164:   DMDAGetBoundingBox(stokes_da,xymin,xymax);

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

1168:   tp_L2     = tu_L2 = tu_H1 = 0.0;
1169:   tint_p_ms = tint_p = 0.0;

1171:   DMDAGetElementsCorners(stokes_da,&sex,&sey,&sez);
1172:   DMDAGetElementsSizes(stokes_da,&mx,&my,&mz);
1173:   for (ek = sez; ek < sez+mz; ek++) {
1174:     for (ej = sey; ej < sey+my; ej++) {
1175:       for (ei = sex; ei < sex+mx; ei++) {
1176:         /* get coords for the element */
1177:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1178:         StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1179:         StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);

1181:         /* evaluate integral */
1182:         for (p = 0; p < ngp; p++) {
1183:           ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1184:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1185:           ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1186:           fac = gp_weight[p]*J_p;

1188:           for (i = 0; i < NODES_PER_EL; i++) {
1189:             tint_p_ms = tint_p_ms+fac*Ni_p[i]*stokes_analytic_e[i].p_dof;
1190:             tint_p    = tint_p   +fac*Ni_p[i]*stokes_e[i].p_dof;
1191:           }
1192:         }
1193:       }
1194:     }
1195:   }
1196:   MPI_Allreduce(&tint_p_ms,&int_p_ms,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1197:   MPI_Allreduce(&tint_p,&int_p,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);

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

1201:   /* remove mine and add manufacture one */
1202:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1203:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);

1205:   {
1206:     PetscInt    k,L,dof;
1207:     PetscScalar *fields;
1208:     DMDAGetInfo(stokes_da,0, 0,0,0, 0,0,0, &dof,0,0,0,0,0);

1210:     VecGetLocalSize(X_local,&L);
1211:     VecGetArray(X_local,&fields);
1212:     for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1213:     VecRestoreArray(X_local,&fields);

1215:     VecGetLocalSize(X,&L);
1216:     VecGetArray(X,&fields);
1217:     for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1218:     VecRestoreArray(X,&fields);
1219:   }

1221:   DMDAVecGetArray(stokes_da,X_local,&stokes);
1222:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1224:   for (ek = sez; ek < sez+mz; ek++) {
1225:     for (ej = sey; ej < sey+my; ej++) {
1226:       for (ei = sex; ei < sex+mx; ei++) {
1227:         /* get coords for the element */
1228:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1229:         StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1230:         StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);

1232:         /* evaluate integral */
1233:         p_e_L2 = 0.0;
1234:         u_e_L2 = 0.0;
1235:         u_e_H1 = 0.0;
1236:         for (p = 0; p < ngp; p++) {
1237:           ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1238:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1239:           ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1240:           fac = gp_weight[p]*J_p;

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

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

1247:             u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1248:             v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1249:             w_error = stokes_e[i].w_dof-stokes_analytic_e[i].w_dof;
1250:             u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error+w_error*w_error);

1252:             u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error              /* du/dx */
1253:                                  +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error
1254:                                  +GNx_p[2][i]*u_error*GNx_p[2][i]*u_error
1255:                                  +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error               /* dv/dx */
1256:                                  +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error
1257:                                  +GNx_p[2][i]*v_error*GNx_p[2][i]*v_error
1258:                                  +GNx_p[0][i]*w_error*GNx_p[0][i]*w_error               /* dw/dx */
1259:                                  +GNx_p[1][i]*w_error*GNx_p[1][i]*w_error
1260:                                  +GNx_p[2][i]*w_error*GNx_p[2][i]*w_error);
1261:           }
1262:         }

1264:         tp_L2 += p_e_L2;
1265:         tu_L2 += u_e_L2;
1266:         tu_H1 += u_e_H1;
1267:       }
1268:     }
1269:   }
1270:   MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1271:   MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1272:   MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1273:   p_L2 = PetscSqrtScalar(p_L2);
1274:   u_L2 = PetscSqrtScalar(u_L2);
1275:   u_H1 = PetscSqrtScalar(u_H1);

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


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

1282:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1283:   VecDestroy(&X_analytic_local);
1284:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1285:   VecDestroy(&X_local);
1286:   return(0);
1287: }

1289: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da,Vec FIELD,const char file_prefix[])
1290: {
1291:   char           vtk_filename[PETSC_MAX_PATH_LEN];
1292:   PetscMPIInt    rank;
1293:   MPI_Comm       comm;
1294:   FILE           *vtk_fp = NULL;
1295:   PetscInt       si,sj,sk,nx,ny,nz,i;
1296:   PetscInt       f,n_fields,N;
1297:   DM             cda;
1298:   Vec            coords;
1299:   Vec            l_FIELD;
1300:   PetscScalar    *_L_FIELD;
1301:   PetscInt       memory_offset;
1302:   PetscScalar    *buffer;


1307:   /* create file name */
1308:   PetscObjectGetComm((PetscObject)da,&comm);
1309:   MPI_Comm_rank(comm,&rank);
1310:   PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"subdomain-%s-p%1.4d.vts",file_prefix,rank);

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

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

1318:   /* coords */
1319:   DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1320:   N    = nx * ny * nz;

1322: #if defined(PETSC_WORDS_BIGENDIAN)
1323:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1324: #else
1325:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1326: #endif
1327:   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);
1328:   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);

1330:   memory_offset = 0;

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

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

1336:   /* copy coordinates */
1337:   DMGetCoordinateDM(da,&cda);
1338:   DMGetCoordinatesLocal(da,&coords);
1339:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"        <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%d\" />\n",memory_offset);
1340:   memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N*3;

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

1344:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      <PointData Scalars=\" ");
1345:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_fields,0,0,0,0,0);
1346:   for (f=0; f<n_fields; f++) {
1347:     const char *field_name;
1348:     DMDAGetFieldName(da,f,&field_name);
1349:     PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"%s ",field_name);
1350:   }
1351:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\">\n");

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

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

1361:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      </PointData>\n");
1362:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </Piece>\n");
1363:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  </StructuredGrid>\n");

1365:   PetscMalloc1(N,&buffer);
1366:   DMGetLocalVector(da,&l_FIELD);
1367:   DMGlobalToLocalBegin(da, FIELD,INSERT_VALUES,l_FIELD);
1368:   DMGlobalToLocalEnd(da,FIELD,INSERT_VALUES,l_FIELD);
1369:   VecGetArray(l_FIELD,&_L_FIELD);

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

1374:   /* write coordinates */
1375:   {
1376:     int         length = sizeof(PetscScalar)*N*3;
1377:     PetscScalar *allcoords;

1379:     fwrite(&length,sizeof(int),1,vtk_fp);
1380:     VecGetArray(coords,&allcoords);
1381:     fwrite(allcoords,sizeof(PetscScalar),3*N,vtk_fp);
1382:     VecRestoreArray(coords,&allcoords);
1383:   }
1384:   /* write fields */
1385:   for (f=0; f<n_fields; f++) {
1386:     int length = sizeof(PetscScalar)*N;
1387:     fwrite(&length,sizeof(int),1,vtk_fp);
1388:     /* load */
1389:     for (i=0; i<N; i++) buffer[i] = _L_FIELD[n_fields*i + f];

1391:     /* write */
1392:     fwrite(buffer,sizeof(PetscScalar),N,vtk_fp);
1393:   }
1394:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\n  </AppendedData>\n");

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

1398:   PetscFree(buffer);
1399:   VecRestoreArray(l_FIELD,&_L_FIELD);
1400:   DMRestoreLocalVector(da,&l_FIELD);

1402:   if (vtk_fp) {
1403:     fclose(vtk_fp);
1404:     vtk_fp = NULL;
1405:   }

1407:   return(0);
1408: }

1410: PetscErrorCode DAViewVTK_write_PieceExtend(FILE *vtk_fp,PetscInt indent_level,DM da,const char local_file_prefix[])
1411: {
1412:   PetscMPIInt    size,rank;
1413:   MPI_Comm       comm;
1414:   const PetscInt *lx,*ly,*lz;
1415:   PetscInt       M,N,P,pM,pN,pP,sum,*olx,*oly,*olz;
1416:   PetscInt       *osx,*osy,*osz,*oex,*oey,*oez;
1417:   PetscInt       i,j,k,II,stencil;

1421:   /* create file name */
1422:   PetscObjectGetComm((PetscObject)da,&comm);
1423:   MPI_Comm_size(comm,&size);
1424:   MPI_Comm_rank(comm,&rank);

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

1429:   /* generate start,end list */
1430:   PetscMalloc1(pM+1,&olx);
1431:   PetscMalloc1(pN+1,&oly);
1432:   PetscMalloc1(pP+1,&olz);
1433:   sum  = 0;
1434:   for (i=0; i<pM; i++) {
1435:     olx[i] = sum;
1436:     sum    = sum + lx[i];
1437:   }
1438:   olx[pM] = sum;
1439:   sum     = 0;
1440:   for (i=0; i<pN; i++) {
1441:     oly[i] = sum;
1442:     sum    = sum + ly[i];
1443:   }
1444:   oly[pN] = sum;
1445:   sum     = 0;
1446:   for (i=0; i<pP; i++) {
1447:     olz[i] = sum;
1448:     sum    = sum + lz[i];
1449:   }
1450:   olz[pP] = sum;

1452:   PetscMalloc1(pM,&osx);
1453:   PetscMalloc1(pN,&osy);
1454:   PetscMalloc1(pP,&osz);
1455:   PetscMalloc1(pM,&oex);
1456:   PetscMalloc1(pN,&oey);
1457:   PetscMalloc1(pP,&oez);
1458:   for (i=0; i<pM; i++) {
1459:     osx[i] = olx[i] - stencil;
1460:     oex[i] = olx[i] + lx[i] + stencil;
1461:     if (osx[i]<0) osx[i]=0;
1462:     if (oex[i]>M) oex[i]=M;
1463:   }

1465:   for (i=0; i<pN; i++) {
1466:     osy[i] = oly[i] - stencil;
1467:     oey[i] = oly[i] + ly[i] + stencil;
1468:     if (osy[i]<0)osy[i]=0;
1469:     if (oey[i]>M)oey[i]=N;
1470:   }
1471:   for (i=0; i<pP; i++) {
1472:     osz[i] = olz[i] - stencil;
1473:     oez[i] = olz[i] + lz[i] + stencil;
1474:     if (osz[i]<0) osz[i]=0;
1475:     if (oez[i]>P) oez[i]=P;
1476:   }

1478:   for (k=0; k<pP; k++) {
1479:     for (j=0; j<pN; j++) {
1480:       for (i=0; i<pM; i++) {
1481:         char     name[PETSC_MAX_PATH_LEN];
1482:         PetscInt procid = i + j*pM + k*pM*pN; /* convert proc(i,j,k) to pid */
1483:         PetscSNPrintf(name,sizeof(name),"subdomain-%s-p%1.4d.vts",local_file_prefix,procid);
1484:         for (II=0; II<indent_level; II++) PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  ");

1486:         PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<Piece Extent=\"%d %d %d %d %d %d\"      Source=\"%s\"/>\n",
1487:                      osx[i],oex[i]-1,
1488:                      osy[j],oey[j]-1,
1489:                      osz[k],oez[k]-1,name);
1490:       }
1491:     }
1492:   }
1493:   PetscFree(olx);
1494:   PetscFree(oly);
1495:   PetscFree(olz);
1496:   PetscFree(osx);
1497:   PetscFree(osy);
1498:   PetscFree(osz);
1499:   PetscFree(oex);
1500:   PetscFree(oey);
1501:   PetscFree(oez);
1502:   return(0);
1503: }

1505: PetscErrorCode DAView_3DVTK_PStructuredGrid(DM da,const char file_prefix[],const char local_file_prefix[])
1506: {
1507:   MPI_Comm       comm;
1508:   PetscMPIInt    size,rank;
1509:   char           vtk_filename[PETSC_MAX_PATH_LEN];
1510:   FILE           *vtk_fp = NULL;
1511:   PetscInt       M,N,P,si,sj,sk,nx,ny,nz;
1512:   PetscInt       i,dofs;

1516:   /* only master generates this file */
1517:   PetscObjectGetComm((PetscObject)da,&comm);
1518:   MPI_Comm_size(comm,&size);
1519:   MPI_Comm_rank(comm,&rank);

1521:   if (rank != 0) return(0);

1523:   /* create file name */
1524:   PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"%s.pvts",file_prefix);
1525:   vtk_fp = fopen(vtk_filename,"w");
1526:   if (!vtk_fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"Cannot open file = %s \n",vtk_filename);

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

1531: #if defined(PETSC_WORDS_BIGENDIAN)
1532:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1533: #else
1534:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1535: #endif

1537:   /* define size of the nodal mesh based on the cell DM */
1538:   DMDAGetInfo(da,0,&M,&N,&P,0,0,0,&dofs,0,0,0,0,0);
1539:   DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1540:   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 */

1542:   /* DUMP THE CELL REFERENCES */
1543:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    <PCellData>\n");
1544:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </PCellData>\n");

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

1550:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    <PPointData>\n");
1551:   for (i=0; i<dofs; i++) {
1552:     const char *fieldname;
1553:     DMDAGetFieldName(da,i,&fieldname);
1554:     PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n",fieldname);
1555:   }
1556:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </PPointData>\n");

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

1561:   /* close the file */
1562:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  </PStructuredGrid>\n");
1563:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");

1565:   if (vtk_fp) {
1566:     fclose(vtk_fp);
1567:     vtk_fp = NULL;
1568:   }
1569:   return(0);
1570: }

1572: PetscErrorCode DAView3DPVTS(DM da, Vec x,const char NAME[])
1573: {
1574:   char           vts_filename[PETSC_MAX_PATH_LEN];
1575:   char           pvts_filename[PETSC_MAX_PATH_LEN];

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

1582:   PetscSNPrintf(pvts_filename,sizeof(pvts_filename),"%s-mesh",NAME);
1583:   DAView_3DVTK_PStructuredGrid(da,pvts_filename,vts_filename);
1584:   return(0);
1585: }

1587: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
1588: {
1590:   PetscReal      norms[4];
1591:   Vec            Br,v,w;
1592:   Mat            A;

1595:   KSPGetOperators(ksp,&A,NULL);
1596:   MatCreateVecs(A,&w,&v);

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

1600:   VecStrideNorm(Br,0,NORM_2,&norms[0]);
1601:   VecStrideNorm(Br,1,NORM_2,&norms[1]);
1602:   VecStrideNorm(Br,2,NORM_2,&norms[2]);
1603:   VecStrideNorm(Br,3,NORM_2,&norms[3]);

1605:   VecDestroy(&v);
1606:   VecDestroy(&w);

1608:   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]);
1609:   return(0);
1610: }

1612: static PetscErrorCode PCMGSetupViaCoarsen(PC pc,DM da_fine)
1613: {
1614:   PetscInt              nlevels,k;
1615:   PETSC_UNUSED PetscInt finest;
1616:   DM                    *da_list,*daclist;
1617:   Mat                   R;
1618:   PetscErrorCode        ierr;

1621:   nlevels = 1;
1622:   PetscOptionsGetInt(NULL,NULL,"-levels",&nlevels,0);

1624:   PetscMalloc1(nlevels,&da_list);
1625:   for (k=0; k<nlevels; k++) da_list[k] = NULL;
1626:   PetscMalloc1(nlevels,&daclist);
1627:   for (k=0; k<nlevels; k++) daclist[k] = NULL;

1629:   /* finest grid is nlevels - 1 */
1630:   finest     = nlevels - 1;
1631:   daclist[0] = da_fine;
1632:   PetscObjectReference((PetscObject)da_fine);
1633:   DMCoarsenHierarchy(da_fine,nlevels-1,&daclist[1]);
1634:   for (k=0; k<nlevels; k++) {
1635:     da_list[k] = daclist[nlevels-1-k];
1636:     DMDASetUniformCoordinates(da_list[k],0.0,1.0,0.0,1.0,0.0,1.0);
1637:   }

1639:   PCMGSetLevels(pc,nlevels,NULL);
1640:   PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
1641:   PCMGSetGalerkin(pc,PC_MG_GALERKIN_PMAT);

1643:   for (k=1; k<nlevels; k++) {
1644:     DMCreateInterpolation(da_list[k-1],da_list[k],&R,NULL);
1645:     PCMGSetInterpolation(pc,k,R);
1646:     MatDestroy(&R);
1647:   }

1649:   /* tidy up */
1650:   for (k=0; k<nlevels; k++) {
1651:     DMDestroy(&da_list[k]);
1652:   }
1653:   PetscFree(da_list);
1654:   PetscFree(daclist);
1655:   return(0);
1656: }

1658: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx,PetscInt my,PetscInt mz)
1659: {
1660:   DM             da_Stokes;
1661:   PetscInt       u_dof,p_dof,dof,stencil_width;
1662:   Mat            A,B;
1663:   DM             vel_cda;
1664:   Vec            vel_coords;
1665:   PetscInt       p;
1666:   Vec            f,X,X1;
1667:   DMDACoor3d     ***_vel_coords;
1668:   PetscInt       its;
1669:   KSP            ksp_S;
1670:   PetscInt       model_definition = 0;
1671:   PetscInt       ei,ej,ek,sex,sey,sez,Imx,Jmy,Kmz;
1672:   CellProperties cell_properties;
1673:   PetscBool      write_output = PETSC_FALSE,resolve= PETSC_FALSE;

1677:   PetscOptionsGetBool(NULL,NULL,"-resolve",&resolve,NULL);
1678:   /* Generate the da for velocity and pressure */
1679:   /* Num nodes in each direction is mx+1, my+1, mz+1 */
1680:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1681:   p_dof         = P_DOFS; /* p - pressure */
1682:   dof           = u_dof+p_dof;
1683:   stencil_width = 1;
1684:   DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1685:                                mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,NULL,&da_Stokes);
1686:   DMSetMatType(da_Stokes,MATAIJ);
1687:   DMSetFromOptions(da_Stokes);
1688:   DMSetUp(da_Stokes);
1689:   DMDASetFieldName(da_Stokes,0,"Vx");
1690:   DMDASetFieldName(da_Stokes,1,"Vy");
1691:   DMDASetFieldName(da_Stokes,2,"Vz");
1692:   DMDASetFieldName(da_Stokes,3,"P");

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

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

1700:   /* interpolate the coordinates to quadrature points */
1701:   DMGetCoordinateDM(da_Stokes,&vel_cda);
1702:   DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1703:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1704:   DMDAGetElementsCorners(da_Stokes,&sex,&sey,&sez);
1705:   DMDAGetElementsSizes(da_Stokes,&Imx,&Jmy,&Kmz);
1706:   for (ek = sez; ek < sez+Kmz; ek++) {
1707:     for (ej = sey; ej < sey+Jmy; ej++) {
1708:       for (ei = sex; ei < sex+Imx; ei++) {
1709:         /* get coords for the element */
1710:         PetscInt               ngp;
1711:         PetscScalar            gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
1712:         PetscScalar            el_coords[NSD*NODES_PER_EL];
1713:         GaussPointCoefficients *cell;

1715:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1716:         GetElementCoords3D(_vel_coords,ei,ej,ek,el_coords);
1717:         ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

1719:         for (p = 0; p < GAUSS_POINTS; p++) {
1720:           PetscScalar xi_p[NSD],Ni_p[NODES_PER_EL];
1721:           PetscScalar gp_x,gp_y,gp_z;
1722:           PetscInt    n;

1724:           xi_p[0] = gp_xi[p][0];
1725:           xi_p[1] = gp_xi[p][1];
1726:           xi_p[2] = gp_xi[p][2];
1727:           ShapeFunctionQ13D_Evaluate(xi_p,Ni_p);

1729:           gp_x = gp_y = gp_z = 0.0;
1730:           for (n = 0; n < NODES_PER_EL; n++) {
1731:             gp_x = gp_x+Ni_p[n]*el_coords[NSD*n];
1732:             gp_y = gp_y+Ni_p[n]*el_coords[NSD*n+1];
1733:             gp_z = gp_z+Ni_p[n]*el_coords[NSD*n+2];
1734:           }
1735:           cell->gp_coords[NSD*p]   = gp_x;
1736:           cell->gp_coords[NSD*p+1] = gp_y;
1737:           cell->gp_coords[NSD*p+2] = gp_z;
1738:         }
1739:       }
1740:     }
1741:   }

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

1745:   switch (model_definition) {
1746:   case 0: /* isoviscous */
1747:     for (ek = sez; ek < sez+Kmz; ek++) {
1748:       for (ej = sey; ej < sey+Jmy; ej++) {
1749:         for (ei = sex; ei < sex+Imx; ei++) {
1750:           GaussPointCoefficients *cell;

1752:           CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1753:           for (p = 0; p < GAUSS_POINTS; p++) {
1754:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1755:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1756:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);

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

1760:             cell->fx[p] = 0.0*coord_x;
1761:             cell->fy[p] = -PetscSinReal(2.2*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1762:             cell->fz[p] = 0.0*coord_z;
1763:             cell->hc[p] = 0.0;
1764:           }
1765:         }
1766:       }
1767:     }
1768:     break;

1770:   case 1: /* manufactured */
1771:     for (ek = sez; ek < sez+Kmz; ek++) {
1772:       for (ej = sey; ej < sey+Jmy; ej++) {
1773:         for (ei = sex; ei < sex+Imx; ei++) {
1774:           PetscReal              eta,Fm[NSD],Fc,pos[NSD];
1775:           GaussPointCoefficients *cell;

1777:           CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1778:           for (p = 0; p < GAUSS_POINTS; p++) {
1779:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1780:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1781:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);

1783:             pos[0] = coord_x;
1784:             pos[1] = coord_y;
1785:             pos[2] = coord_z;

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

1790:             cell->fx[p] = Fm[0];
1791:             cell->fy[p] = Fm[1];
1792:             cell->fz[p] = Fm[2];
1793:             cell->hc[p] = Fc;
1794:           }
1795:         }
1796:       }
1797:     }
1798:     break;

1800:   case 2: /* solcx */
1801:     for (ek = sez; ek < sez+Kmz; ek++) {
1802:       for (ej = sey; ej < sey+Jmy; ej++) {
1803:         for (ei = sex; ei < sex+Imx; ei++) {
1804:           GaussPointCoefficients *cell;

1806:           CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1807:           for (p = 0; p < GAUSS_POINTS; p++) {
1808:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1809:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1810:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);

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

1814:             cell->fx[p] = 0.0;
1815:             cell->fy[p] = -PetscSinReal(3.0*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1816:             cell->fz[p] = 0.0*coord_z;
1817:             cell->hc[p] = 0.0;
1818:           }
1819:         }
1820:       }
1821:     }
1822:     break;

1824:   case 3: /* sinker */
1825:     for (ek = sez; ek < sez+Kmz; ek++) {
1826:       for (ej = sey; ej < sey+Jmy; ej++) {
1827:         for (ei = sex; ei < sex+Imx; ei++) {
1828:           GaussPointCoefficients *cell;

1830:           CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1831:           for (p = 0; p < GAUSS_POINTS; p++) {
1832:             PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p]);
1833:             PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);
1834:             PetscReal zp = PetscRealPart(cell->gp_coords[NSD*p+2]);

1836:             cell->eta[p] = 1.0e-2;
1837:             cell->fx[p]  = 0.0;
1838:             cell->fy[p]  = 0.0;
1839:             cell->fz[p]  = 0.0;
1840:             cell->hc[p]  = 0.0;

1842:             if ((PetscAbs(xp-0.5) < 0.2) && (PetscAbs(yp-0.5) < 0.2) && (PetscAbs(zp-0.5) < 0.2)) {
1843:               cell->eta[p] = 1.0;
1844:               cell->fz[p]  = 1.0;
1845:             }

1847:           }
1848:         }
1849:       }
1850:     }
1851:     break;

1853:   case 4: /* subdomain jumps */
1854:     for (ek = sez; ek < sez+Kmz; ek++) {
1855:       for (ej = sey; ej < sey+Jmy; ej++) {
1856:         for (ei = sex; ei < sex+Imx; ei++) {
1857:           PetscReal              opts_mag,opts_eta0;
1858:           PetscInt               px,py,pz;
1859:           PetscBool              jump;
1860:           PetscMPIInt            rr;
1861:           GaussPointCoefficients *cell;

1863:           opts_mag  = 1.0;
1864:           opts_eta0 = 1.e-2;

1866:           PetscOptionsGetReal(NULL,NULL,"-jump_eta0",&opts_eta0,NULL);
1867:           PetscOptionsGetReal(NULL,NULL,"-jump_magnitude",&opts_mag,NULL);
1868:           DMDAGetInfo(da_Stokes,NULL,NULL,NULL,NULL,&px,&py,&pz,NULL,NULL,NULL,NULL,NULL,NULL);
1869:           rr   = PetscGlobalRank%(px*py);
1870:           if (px%2) jump = (PetscBool)(rr%2);
1871:           else jump = (PetscBool)((rr/px)%2 ? rr%2 : !(rr%2));
1872:           rr = PetscGlobalRank/(px*py);
1873:           if (rr%2) jump = (PetscBool)!jump;
1874:           CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1875:           for (p = 0; p < GAUSS_POINTS; p++) {
1876:             PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p]);
1877:             PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);

1879:             cell->eta[p] = jump ? PetscPowReal(10.0,opts_mag) : opts_eta0;
1880:             cell->fx[p]  = 0.0;
1881:             cell->fy[p]  = -PetscSinReal(2.2*PETSC_PI*yp)*PetscCosReal(1.0*PETSC_PI*xp);
1882:             cell->fz[p]  = 0.0;
1883:             cell->hc[p]  = 0.0;
1884:           }
1885:         }
1886:       }
1887:     }
1888:     break;
1889:   default:
1890:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"No default model is supported. Choose either -model {0,1,2,3}");
1891:     break;
1892:   }

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

1896:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1897:   DMCreateMatrix(da_Stokes,&A);
1898:   DMCreateMatrix(da_Stokes,&B);
1899:   DMCreateGlobalVector(da_Stokes,&X);
1900:   DMCreateGlobalVector(da_Stokes,&f);

1902:   /* assemble A11 */
1903:   MatZeroEntries(A);
1904:   MatZeroEntries(B);
1905:   VecZeroEntries(f);

1907:   AssembleA_Stokes(A,da_Stokes,cell_properties);
1908:   MatViewFromOptions(A,NULL,"-amat_view");
1909:   AssembleA_PCStokes(B,da_Stokes,cell_properties);
1910:   MatViewFromOptions(B,NULL,"-bmat_view");
1911:   /* build force vector */
1912:   AssembleF_Stokes(f,da_Stokes,cell_properties);

1914:   /* SOLVE */
1915:   KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1916:   KSPSetOptionsPrefix(ksp_S,"stokes_"); /* stokes */
1917:   KSPSetOperators(ksp_S,A,B);
1918:   KSPSetFromOptions(ksp_S);

1920:   {
1921:     PC             pc;
1922:     const PetscInt ufields[] = {0,1,2},pfields[] = {3};
1923:     KSPGetPC(ksp_S,&pc);
1924:     PCFieldSplitSetBlockSize(pc,4);
1925:     PCFieldSplitSetFields(pc,"u",3,ufields,ufields);
1926:     PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1927:   }

1929:   {
1930:     PC        pc;
1931:     PetscBool same = PETSC_FALSE;
1932:     KSPGetPC(ksp_S,&pc);
1933:     PetscObjectTypeCompare((PetscObject)pc,PCMG,&same);
1934:     if (same) {
1935:       PCMGSetupViaCoarsen(pc,da_Stokes);
1936:     }
1937:   }

1939:   {
1940:     PC        pc;
1941:     PetscBool same = PETSC_FALSE;
1942:     KSPGetPC(ksp_S,&pc);
1943:     PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&same);
1944:     if (same) {
1945:       KSPSetOperators(ksp_S,A,A);
1946:     }
1947:   }

1949:   {
1950:     PetscBool stokes_monitor = PETSC_FALSE;
1951:     PetscOptionsGetBool(NULL,NULL,"-stokes_ksp_monitor_blocks",&stokes_monitor,0);
1952:     if (stokes_monitor) {
1953:       KSPMonitorSet(ksp_S,KSPMonitorStokesBlocks,NULL,NULL);
1954:     }
1955:   }

1957:   if (resolve) {
1958:     /* Test changing matrix data structure and resolve */
1959:     VecDuplicate(X,&X1);
1960:   }

1962:   KSPSolve(ksp_S,f,X);
1963:   if (resolve) {
1964:     Mat C;
1965:     MatDuplicate(A,MAT_COPY_VALUES,&C);
1966:     KSPSetOperators(ksp_S,C,C);
1967:     KSPSolve(ksp_S,f,X1);
1968:     MatDestroy(&C);
1969:     VecDestroy(&X1);
1970:   }

1972:   PetscOptionsGetBool(NULL,NULL,"-write_pvts",&write_output,NULL);
1973:   if (write_output) {DAView3DPVTS(da_Stokes,X,"up");}
1974:   {
1975:     PetscBool flg = PETSC_FALSE;
1976:     char      filename[PETSC_MAX_PATH_LEN];
1977:     PetscOptionsGetString(NULL,NULL,"-write_binary",filename,sizeof(filename),&flg);
1978:     if (flg) {
1979:       PetscViewer viewer;
1980:       /* PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer); */
1981:       PetscViewerVTKOpen(PETSC_COMM_WORLD,"ex42.vts",FILE_MODE_WRITE,&viewer);
1982:       VecView(X,viewer);
1983:       PetscViewerDestroy(&viewer);
1984:     }
1985:   }
1986:   KSPGetIterationNumber(ksp_S,&its);

1988:   /* verify */
1989:   if (model_definition == 1) {
1990:     DM  da_Stokes_analytic;
1991:     Vec X_analytic;

1993:     DMDACreateManufacturedSolution(mx,my,mz,&da_Stokes_analytic,&X_analytic);
1994:     if (write_output) {
1995:       DAView3DPVTS(da_Stokes_analytic,X_analytic,"ms");
1996:     }
1997:     DMDAIntegrateErrors3D(da_Stokes_analytic,X,X_analytic);
1998:     if (write_output) {
1999:       DAView3DPVTS(da_Stokes,X,"up2");
2000:     }
2001:     DMDestroy(&da_Stokes_analytic);
2002:     VecDestroy(&X_analytic);
2003:   }

2005:   KSPDestroy(&ksp_S);
2006:   VecDestroy(&X);
2007:   VecDestroy(&f);
2008:   MatDestroy(&A);
2009:   MatDestroy(&B);

2011:   CellPropertiesDestroy(&cell_properties);
2012:   DMDestroy(&da_Stokes);
2013:   return(0);
2014: }

2016: int main(int argc,char **args)
2017: {
2019:   PetscInt       mx,my,mz;

2021:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
2022:   mx   = my = mz = 10;
2023:   PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
2024:   my   = mx; mz = mx;
2025:   PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
2026:   PetscOptionsGetInt(NULL,NULL,"-mz",&mz,NULL);
2027:   solve_stokes_3d_coupled(mx,my,mz);
2028:   PetscFinalize();
2029:   return ierr;
2030: }


2033: /*TEST

2035:    test:
2036:       args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type lu

2038:    test:
2039:       suffix: 2
2040:       nsize: 3
2041:       args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type redundant -stokes_redundant_pc_type lu

2043:    test:
2044:       suffix: bddc_stokes
2045:       nsize: 6
2046:       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

2048:    test:
2049:       suffix: bddc_stokes_deluxe
2050:       nsize: 6
2051:       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

2053:    test:
2054:       requires: !single
2055:       suffix: bddc_stokes_subdomainjump_deluxe
2056:       nsize: 9
2057:       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

2059:    test:
2060:       requires: !single
2061:       suffix: 3
2062:       args: -stokes_ksp_converged_reason -stokes_pc_type fieldsplit -resolve

2064:    test:
2065:       suffix: tut_1
2066:       nsize: 4
2067:       requires: !single
2068:       args: -stokes_ksp_monitor

2070:    test:
2071:       suffix: tut_2
2072:       nsize: 4
2073:       requires: !single
2074:       args: -stokes_ksp_monitor -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_type schur

2076:    test:
2077:       suffix: tut_3
2078:       nsize: 4
2079:       requires: !single
2080:       args: -mx 20 -stokes_ksp_monitor -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_type schur

2082: TEST*/