Actual source code: ex42.c

petsc-master 2017-07-21
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) \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: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz);

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

 59:   PetscNew(&cells);

 61:   DMDAGetElementCorners(da_stokes,&sex,&sey,&sez,&mx,&my,&mz);

 63:   cells->mx     = mx;
 64:   cells->my     = my;
 65:   cells->mz     = mz;
 66:   cells->ncells = mx * my * mz;
 67:   cells->sex    = sex;
 68:   cells->sey    = sey;
 69:   cells->sez    = sez;

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

 73:   *C = cells;
 74:   return(0);
 75: }

 77: PetscErrorCode CellPropertiesDestroy(CellProperties *C)
 78: {
 80:   CellProperties cells;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

215:   matrix_inverse_3x3(JJ,iJ);

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

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

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

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

239:   gp_weight[0] = 1.0;
240:   gp_weight[1] = 1.0;
241:   gp_weight[2] = 1.0;
242:   gp_weight[3] = 1.0;

244:   gp_weight[4] = 1.0;
245:   gp_weight[5] = 1.0;
246:   gp_weight[6] = 1.0;
247:   gp_weight[7] = 1.0;
248: }

250: /* procs to the left claim the ghost node as their element */
251: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
252: {
253:   PetscInt       m,n,p,M,N,P;
254:   PetscInt       sx,sy,sz;

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

261:   if (mxl) {
262:     *mxl = m;
263:     if ((sx+m) == M) *mxl = m-1;  /* last proc */
264:   }
265:   if (myl) {
266:     *myl = n;
267:     if ((sy+n) == N) *myl = n-1;  /* last proc */
268:   }
269:   if (mzl) {
270:     *mzl = p;
271:     if ((sz+p) == P) *mzl = p-1;  /* last proc */
272:   }
273:   return(0);
274: }

276: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
277: {
278:   PetscInt       si,sj,sk;

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

284:   if (sx) {
285:     *sx = si;
286:     if (si != 0) *sx = si+1;
287:   }
288:   if (sy) {
289:     *sy = sj;
290:     if (sj != 0) *sy = sj+1;
291:   }
292:   if (sz) {
293:     *sz = sk;
294:     if (sk != 0) *sz = sk+1;
295:   }
296:   DMDAGetLocalElementSize(da,mx,my,mz);
297:   return(0);
298: }

300: /*
301:  i,j are the element indices
302:  The unknown is a vector quantity.
303:  The s[].c is used to indicate the degree of freedom.
304:  */
305: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j,PetscInt k)
306: {
307:   PetscInt n;

310:   /* velocity */
311:   n = 0;
312:   /* node 0 */
313:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++; /* Vx0 */
314:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++; /* Vy0 */
315:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++; /* Vz0 */

317:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
318:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
319:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;

321:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
322:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
323:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;

325:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++;
326:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++;
327:   s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++;

329:   /* */
330:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++; /* Vx4 */
331:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++; /* Vy4 */
332:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++; /* Vz4 */

334:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
335:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
336:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;

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

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

346:   /* pressure */
347:   n = 0;

349:   s_p[n].i = i;   s_p[n].j = j;   s_p[n].k = k; s_p[n].c = 3; n++; /* P0 */
350:   s_p[n].i = i;   s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
351:   s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
352:   s_p[n].i = i+1; s_p[n].j = j;   s_p[n].k = k; s_p[n].c = 3; n++;

354:   s_p[n].i = i;   s_p[n].j = j;   s_p[n].k = k+1; s_p[n].c = 3; n++; /* P0 */
355:   s_p[n].i = i;   s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
356:   s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
357:   s_p[n].i = i+1; s_p[n].j = j;   s_p[n].k = k+1; s_p[n].c = 3; n++;
358:   return(0);
359: }

361: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords,PetscInt i,PetscInt j,PetscInt k,PetscScalar el_coord[])
362: {
364:   /* get coords for the element */
365:   el_coord[0] = coords[k][j][i].x;
366:   el_coord[1] = coords[k][j][i].y;
367:   el_coord[2] = coords[k][j][i].z;

369:   el_coord[3] = coords[k][j+1][i].x;
370:   el_coord[4] = coords[k][j+1][i].y;
371:   el_coord[5] = coords[k][j+1][i].z;

373:   el_coord[6] = coords[k][j+1][i+1].x;
374:   el_coord[7] = coords[k][j+1][i+1].y;
375:   el_coord[8] = coords[k][j+1][i+1].z;

377:   el_coord[9]  = coords[k][j][i+1].x;
378:   el_coord[10] = coords[k][j][i+1].y;
379:   el_coord[11] = coords[k][j][i+1].z;

381:   el_coord[12] = coords[k+1][j][i].x;
382:   el_coord[13] = coords[k+1][j][i].y;
383:   el_coord[14] = coords[k+1][j][i].z;

385:   el_coord[15] = coords[k+1][j+1][i].x;
386:   el_coord[16] = coords[k+1][j+1][i].y;
387:   el_coord[17] = coords[k+1][j+1][i].z;

389:   el_coord[18] = coords[k+1][j+1][i+1].x;
390:   el_coord[19] = coords[k+1][j+1][i+1].y;
391:   el_coord[20] = coords[k+1][j+1][i+1].z;

393:   el_coord[21] = coords[k+1][j][i+1].x;
394:   el_coord[22] = coords[k+1][j][i+1].y;
395:   el_coord[23] = coords[k+1][j][i+1].z;
396:   return(0);
397: }

399: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field,PetscInt i,PetscInt j,PetscInt k,StokesDOF nodal_fields[])
400: {
402:   /* get the nodal fields for u */
403:   nodal_fields[0].u_dof = field[k][j][i].u_dof;
404:   nodal_fields[0].v_dof = field[k][j][i].v_dof;
405:   nodal_fields[0].w_dof = field[k][j][i].w_dof;

407:   nodal_fields[1].u_dof = field[k][j+1][i].u_dof;
408:   nodal_fields[1].v_dof = field[k][j+1][i].v_dof;
409:   nodal_fields[1].w_dof = field[k][j+1][i].w_dof;

411:   nodal_fields[2].u_dof = field[k][j+1][i+1].u_dof;
412:   nodal_fields[2].v_dof = field[k][j+1][i+1].v_dof;
413:   nodal_fields[2].w_dof = field[k][j+1][i+1].w_dof;

415:   nodal_fields[3].u_dof = field[k][j][i+1].u_dof;
416:   nodal_fields[3].v_dof = field[k][j][i+1].v_dof;
417:   nodal_fields[3].w_dof = field[k][j][i+1].w_dof;

419:   nodal_fields[4].u_dof = field[k+1][j][i].u_dof;
420:   nodal_fields[4].v_dof = field[k+1][j][i].v_dof;
421:   nodal_fields[4].w_dof = field[k+1][j][i].w_dof;

423:   nodal_fields[5].u_dof = field[k+1][j+1][i].u_dof;
424:   nodal_fields[5].v_dof = field[k+1][j+1][i].v_dof;
425:   nodal_fields[5].w_dof = field[k+1][j+1][i].w_dof;

427:   nodal_fields[6].u_dof = field[k+1][j+1][i+1].u_dof;
428:   nodal_fields[6].v_dof = field[k+1][j+1][i+1].v_dof;
429:   nodal_fields[6].w_dof = field[k+1][j+1][i+1].w_dof;

431:   nodal_fields[7].u_dof = field[k+1][j][i+1].u_dof;
432:   nodal_fields[7].v_dof = field[k+1][j][i+1].v_dof;
433:   nodal_fields[7].w_dof = field[k+1][j][i+1].w_dof;

435:   /* get the nodal fields for p */
436:   nodal_fields[0].p_dof = field[k][j][i].p_dof;
437:   nodal_fields[1].p_dof = field[k][j+1][i].p_dof;
438:   nodal_fields[2].p_dof = field[k][j+1][i+1].p_dof;
439:   nodal_fields[3].p_dof = field[k][j][i+1].p_dof;

441:   nodal_fields[4].p_dof = field[k+1][j][i].p_dof;
442:   nodal_fields[5].p_dof = field[k+1][j+1][i].p_dof;
443:   nodal_fields[6].p_dof = field[k+1][j+1][i+1].p_dof;
444:   nodal_fields[7].p_dof = field[k+1][j][i+1].p_dof;
445:   return(0);
446: }

448: 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)
449: {
450:   PetscInt              ij;
451:   PETSC_UNUSED PetscInt r,c,nr,nc;

453:   nr = w_NPE*w_dof;
454:   nc = u_NPE*u_dof;

456:   r = w_dof*wi+wd;
457:   c = u_dof*ui+ud;

459:   ij = r*nc+c;

461:   return ij;
462: }

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

469:   for (n = 0; n<NODES_PER_EL; n++) {
470:     II = u_eqn[NSD*n].i;
471:     J = u_eqn[NSD*n].j;
472:     K = u_eqn[NSD*n].k;

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

476:     II = u_eqn[NSD*n+1].i;
477:     J = u_eqn[NSD*n+1].j;
478:     K = u_eqn[NSD*n+1].k;

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

482:     II = u_eqn[NSD*n+2].i;
483:     J = u_eqn[NSD*n+2].j;
484:     K = u_eqn[NSD*n+2].k;
485:     fields_F[K][J][II].w_dof = fields_F[K][J][II].w_dof+Fe_u[NSD*n+2];

487:     II = p_eqn[n].i;
488:     J = p_eqn[n].j;
489:     K = p_eqn[n].k;

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

493:   }
494:   return(0);
495: }

497: static void FormStressOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
498: {
499:   PetscInt       ngp;
500:   PetscScalar    gp_xi[GAUSS_POINTS][NSD];
501:   PetscScalar    gp_weight[GAUSS_POINTS];
502:   PetscInt       p,i,j,k;
503:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
504:   PetscScalar    J_p,tildeD[6];
505:   PetscScalar    B[6][U_DOFS*NODES_PER_EL];
506:   const PetscInt nvdof = U_DOFS*NODES_PER_EL;

508:   /* define quadrature rule */
509:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

511:   /* evaluate integral */
512:   for (p = 0; p < ngp; p++) {
513:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
514:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);

516:     for (i = 0; i < NODES_PER_EL; i++) {
517:       PetscScalar d_dx_i = GNx_p[0][i];
518:       PetscScalar d_dy_i = GNx_p[1][i];
519:       PetscScalar d_dz_i = GNx_p[2][i];

521:       B[0][3*i] = d_dx_i; B[0][3*i+1] = 0.0;     B[0][3*i+2] = 0.0;
522:       B[1][3*i] = 0.0;    B[1][3*i+1] = d_dy_i;  B[1][3*i+2] = 0.0;
523:       B[2][3*i] = 0.0;    B[2][3*i+1] = 0.0;     B[2][3*i+2] = d_dz_i;

525:       B[3][3*i] = d_dy_i; B[3][3*i+1] = d_dx_i;  B[3][3*i+2] = 0.0;   /* e_xy */
526:       B[4][3*i] = d_dz_i; B[4][3*i+1] = 0.0;     B[4][3*i+2] = d_dx_i; /* e_xz */
527:       B[5][3*i] = 0.0;    B[5][3*i+1] = d_dz_i;  B[5][3*i+2] = d_dy_i; /* e_yz */
528:     }


531:     tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
532:     tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
533:     tildeD[2] = 2.0*gp_weight[p]*J_p*eta[p];

535:     tildeD[3] =     gp_weight[p]*J_p*eta[p];
536:     tildeD[4] =     gp_weight[p]*J_p*eta[p];
537:     tildeD[5] =     gp_weight[p]*J_p*eta[p];

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

553:   }
554:   /* fill lower triangular part */
555: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
556:   for (i = 0; i < nvdof; i++) {
557:     for (j = i; j < nvdof; j++) {
558:       Ke[j*nvdof+i] = Ke[i*nvdof+j];
559:     }
560:   }
561: #endif
562: }

564: static void FormGradientOperatorQ13D(PetscScalar Ke[],PetscScalar coords[])
565: {
566:   PetscInt    ngp;
567:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
568:   PetscScalar gp_weight[GAUSS_POINTS];
569:   PetscInt    p,i,j,di;
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;

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;

584:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
585:       for (di = 0; di < NSD; di++) { /* u dofs */
586:         for (j = 0; j < NODES_PER_EL; j++) {  /* p nodes, p dofs = 1 (ie no loop) */
587:           PetscInt IJ;
588:           IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,3,j,0,NODES_PER_EL,1);

590:           Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
591:         }
592:       }
593:     }
594:   }
595: }

597: static void FormDivergenceOperatorQ13D(PetscScalar De[],PetscScalar coords[])
598: {
599:   PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
600:   PetscInt    i,j;
601:   PetscInt    nr_g,nc_g;

603:   PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
604:   FormGradientOperatorQ13D(Ge,coords);

606:   nr_g = U_DOFS*NODES_PER_EL;
607:   nc_g = P_DOFS*NODES_PER_EL;

609:   for (i = 0; i < nr_g; i++) {
610:     for (j = 0; j < nc_g; j++) {
611:       De[nr_g*j+i] = Ge[nc_g*i+j];
612:     }
613:   }
614: }

616: static void FormStabilisationOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
617: {
618:   PetscInt    ngp;
619:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
620:   PetscScalar gp_weight[GAUSS_POINTS];
621:   PetscInt    p,i,j;
622:   PetscScalar Ni_p[NODES_PER_EL];
623:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
624:   PetscScalar J_p,fac,eta_avg;

626:   /* define quadrature rule */
627:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

629:   /* evaluate integral */
630:   for (p = 0; p < ngp; p++) {
631:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
632:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
633:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
634:     fac = gp_weight[p]*J_p;
635:     /*
636:      for (i = 0; i < NODES_PER_EL; i++) {
637:      for (j = i; j < NODES_PER_EL; j++) {
638:      Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
639:      }
640:      }
641:      */

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

650:   /* scale */
651:   eta_avg = 0.0;
652:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
653:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
654:   fac     = 1.0/eta_avg;

656:   /*
657:    for (i = 0; i < NODES_PER_EL; i++) {
658:    for (j = i; j < NODES_PER_EL; j++) {
659:    Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
660:    #if defined(ASSEMBLE_LOWER_TRIANGULAR)
661:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
662:    #endif
663:    }
664:    }
665:    */

667:   for (i = 0; i < NODES_PER_EL; i++) {
668:     for (j = 0; j < NODES_PER_EL; j++) {
669:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
670:     }
671:   }
672: }

674: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
675: {
676:   PetscInt    ngp;
677:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
678:   PetscScalar gp_weight[GAUSS_POINTS];
679:   PetscInt    p,i,j;
680:   PetscScalar Ni_p[NODES_PER_EL];
681:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
682:   PetscScalar J_p,fac,eta_avg;

684:   /* define quadrature rule */
685:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

687:   /* evaluate integral */
688:   for (p = 0; p < ngp; p++) {
689:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
690:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
691:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
692:     fac = gp_weight[p]*J_p;

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

702:     for (i = 0; i < NODES_PER_EL; i++) {
703:       for (j = 0; j < NODES_PER_EL; j++) {
704:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
705:       }
706:     }
707:   }

709:   /* scale */
710:   eta_avg = 0.0;
711:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
712:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
713:   fac     = 1.0/eta_avg;
714:   /*
715:    for (i = 0; i < NODES_PER_EL; i++) {
716:    for (j = i; j < NODES_PER_EL; j++) {
717:    Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
718:    #if defined(ASSEMBLE_LOWER_TRIANGULAR)
719:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
720:    #endif
721:    }
722:    }
723:    */

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

732: static void FormMomentumRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[],PetscScalar fz[])
733: {
734:   PetscInt    ngp;
735:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
736:   PetscScalar gp_weight[GAUSS_POINTS];
737:   PetscInt    p,i;
738:   PetscScalar Ni_p[NODES_PER_EL];
739:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
740:   PetscScalar J_p,fac;

742:   /* define quadrature rule */
743:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

745:   /* evaluate integral */
746:   for (p = 0; p < ngp; p++) {
747:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
748:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
749:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
750:     fac = gp_weight[p]*J_p;

752:     for (i = 0; i < NODES_PER_EL; i++) {
753:       Fe[NSD*i]   -= fac*Ni_p[i]*fx[p];
754:       Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
755:       Fe[NSD*i+2] -= fac*Ni_p[i]*fz[p];
756:     }
757:   }
758: }

760: static void FormContinuityRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar hc[])
761: {
762:   PetscInt    ngp;
763:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
764:   PetscScalar gp_weight[GAUSS_POINTS];
765:   PetscInt    p,i;
766:   PetscScalar Ni_p[NODES_PER_EL];
767:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
768:   PetscScalar J_p,fac;

770:   /* define quadrature rule */
771:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

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

780:     for (i = 0; i < NODES_PER_EL; i++) Fe[i] -= fac*Ni_p[i]*hc[p];
781:   }
782: }

784: #define _ZERO_ROWCOL_i(A,i) {                   \
785:     PetscInt    KK;                             \
786:     PetscScalar tmp = A[24*(i)+(i)];            \
787:     for (KK=0;KK<24;KK++) A[24*(i)+KK]=0.0;     \
788:     for (KK=0;KK<24;KK++) A[24*KK+(i)]=0.0;     \
789:     A[24*(i)+(i)] = tmp;}                       \

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

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

799: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,CellProperties cell_properties)
800: {
801:   DM                     cda;
802:   Vec                    coords;
803:   DMDACoor3d             ***_coords;
804:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS];
805:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS];
806:   PetscInt               sex,sey,sez,mx,my,mz;
807:   PetscInt               ei,ej,ek;
808:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
809:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
810:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
811:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
812:   PetscScalar            el_coords[NODES_PER_EL*NSD];
813:   GaussPointCoefficients *props;
814:   PetscScalar            *prop_eta;
815:   PetscInt               n,M,N,P;
816:   PetscErrorCode         ierr;

819:   DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
820:   /* setup for coords */
821:   DMGetCoordinateDM(stokes_da,&cda);
822:   DMGetCoordinatesLocal(stokes_da,&coords);
823:   DMDAVecGetArray(cda,coords,&_coords);

825:   DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
826:   for (ek = sez; ek < sez+mz; ek++) {
827:     for (ej = sey; ej < sey+my; ej++) {
828:       for (ei = sex; ei < sex+mx; ei++) {
829:         /* get coords for the element */
830:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
831:         /* get cell properties */
832:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
833:         /* get coefficients for the element */
834:         prop_eta = props->eta;

836:         /* initialise element stiffness matrix */
837:         PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
838:         PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
839:         PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
840:         PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);

842:         /* form element stiffness matrix */
843:         FormStressOperatorQ13D(Ae,el_coords,prop_eta);
844:         FormGradientOperatorQ13D(Ge,el_coords);
845:         /*#if defined(ASSEMBLE_LOWER_TRIANGULAR)*/
846:         FormDivergenceOperatorQ13D(De,el_coords);
847:         /*#endif*/
848:         FormStabilisationOperatorQ13D(Ce,el_coords,prop_eta);

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

853:         for (n=0; n<NODES_PER_EL; n++) {
854:           if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
855:             _ZERO_ROWCOL_i(Ae,3*n);
856:             _ZERO_ROW_i   (Ge,3*n);
857:             _ZERO_COL_i   (De,3*n);
858:           }

860:           if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
861:             _ZERO_ROWCOL_i(Ae,3*n+1);
862:             _ZERO_ROW_i   (Ge,3*n+1);
863:             _ZERO_COL_i   (De,3*n+1);
864:           }

866:           if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
867:             _ZERO_ROWCOL_i(Ae,3*n+2);
868:             _ZERO_ROW_i   (Ge,3*n+2);
869:             _ZERO_COL_i   (De,3*n+2);
870:           }
871:         }
872:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
873:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
874:         MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
875:         MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
876:       }
877:     }
878:   }
879:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
880:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

884:   return(0);
885: }

887: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,CellProperties cell_properties)
888: {
889:   DM                     cda;
890:   Vec                    coords;
891:   DMDACoor3d             ***_coords;
892:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS];
893:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS];
894:   PetscInt               sex,sey,sez,mx,my,mz;
895:   PetscInt               ei,ej,ek;
896:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
897:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
898:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
899:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
900:   PetscScalar            el_coords[NODES_PER_EL*NSD];
901:   GaussPointCoefficients *props;
902:   PetscScalar            *prop_eta;
903:   PetscInt               n,M,N,P;
904:   PetscErrorCode         ierr;

907:   DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
908:   /* setup for coords */
909:   DMGetCoordinateDM(stokes_da,&cda);
910:   DMGetCoordinatesLocal(stokes_da,&coords);
911:   DMDAVecGetArray(cda,coords,&_coords);

913:   DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
914:   for (ek = sez; ek < sez+mz; ek++) {
915:     for (ej = sey; ej < sey+my; ej++) {
916:       for (ei = sex; ei < sex+mx; ei++) {
917:         /* get coords for the element */
918:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
919:         /* get cell properties */
920:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
921:         /* get coefficients for the element */
922:         prop_eta = props->eta;

924:         /* initialise element stiffness matrix */
925:         PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
926:         PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
927:         PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
928:         PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);

930:         /* form element stiffness matrix */
931:         FormStressOperatorQ13D(Ae,el_coords,prop_eta);
932:         FormGradientOperatorQ13D(Ge,el_coords);
933:         /* FormDivergenceOperatorQ13D(De,el_coords); */
934:         FormScaledMassMatrixOperatorQ13D(Ce,el_coords,prop_eta);

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

939:         for (n=0; n<NODES_PER_EL; n++) {
940:           if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
941:             _ZERO_ROWCOL_i(Ae,3*n);
942:             _ZERO_ROW_i   (Ge,3*n);
943:             _ZERO_COL_i   (De,3*n);
944:           }

946:           if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
947:             _ZERO_ROWCOL_i(Ae,3*n+1);
948:             _ZERO_ROW_i   (Ge,3*n+1);
949:             _ZERO_COL_i   (De,3*n+1);
950:           }

952:           if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
953:             _ZERO_ROWCOL_i(Ae,3*n+2);
954:             _ZERO_ROW_i   (Ge,3*n+2);
955:             _ZERO_COL_i   (De,3*n+2);
956:           }
957:         }
958:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
959:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
960:         /*MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);*/
961:         MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
962:       }
963:     }
964:   }
965:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
966:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

968:   DMDAVecRestoreArray(cda,coords,&_coords);
969:   return(0);
970: }

972: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,CellProperties cell_properties)
973: {
974:   DM                     cda;
975:   Vec                    coords;
976:   DMDACoor3d             ***_coords;
977:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS];
978:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS];
979:   PetscInt               sex,sey,sez,mx,my,mz;
980:   PetscInt               ei,ej,ek;
981:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
982:   PetscScalar            He[NODES_PER_EL*P_DOFS];
983:   PetscScalar            el_coords[NODES_PER_EL*NSD];
984:   GaussPointCoefficients *props;
985:   PetscScalar            *prop_fx,*prop_fy,*prop_fz,*prop_hc;
986:   Vec                    local_F;
987:   StokesDOF              ***ff;
988:   PetscInt               n,M,N,P;
989:   PetscErrorCode         ierr;

992:   DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
993:   /* setup for coords */
994:   DMGetCoordinateDM(stokes_da,&cda);
995:   DMGetCoordinatesLocal(stokes_da,&coords);
996:   DMDAVecGetArray(cda,coords,&_coords);

998:   /* get acces to the vector */
999:   DMGetLocalVector(stokes_da,&local_F);
1000:   VecZeroEntries(local_F);
1001:   DMDAVecGetArray(stokes_da,local_F,&ff);

1003:   DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
1004:   for (ek = sez; ek < sez+mz; ek++) {
1005:     for (ej = sey; ej < sey+my; ej++) {
1006:       for (ei = sex; ei < sex+mx; ei++) {
1007:         /* get coords for the element */
1008:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1009:         /* get cell properties */
1010:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
1011:         /* get coefficients for the element */
1012:         prop_fx = props->fx;
1013:         prop_fy = props->fy;
1014:         prop_fz = props->fz;
1015:         prop_hc = props->hc;

1017:         /* initialise element stiffness matrix */
1018:         PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
1019:         PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);

1021:         /* form element stiffness matrix */
1022:         FormMomentumRhsQ13D(Fe,el_coords,prop_fx,prop_fy,prop_fz);
1023:         FormContinuityRhsQ13D(He,el_coords,prop_hc);

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

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

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

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

1036:         DMDASetValuesLocalStencil3D_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
1037:       }
1038:     }
1039:   }
1040:   DMDAVecRestoreArray(stokes_da,local_F,&ff);
1041:   DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
1042:   DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
1043:   DMRestoreLocalVector(stokes_da,&local_F);

1045:   DMDAVecRestoreArray(cda,coords,&_coords);
1046:   return(0);
1047: }

1049: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta,PetscReal *MX,PetscReal *MY,PetscReal *MZ)
1050: {
1051:   *theta = 0.0;
1052:   *MX    = 2.0 * PETSC_PI;
1053:   *MY    = 2.0 * PETSC_PI;
1054:   *MZ    = 2.0 * PETSC_PI;
1055: }
1056: static void evaluate_MS_FrankKamentski(PetscReal pos[],PetscReal v[],PetscReal *p,PetscReal *eta,PetscReal Fm[],PetscReal *Fc)
1057: {
1058:   PetscReal x,y,z;
1059:   PetscReal theta,MX,MY,MZ;

1061:   evaluate_MS_FrankKamentski_constants(&theta,&MX,&MY,&MZ);
1062:   x = pos[0];
1063:   y = pos[1];
1064:   z = pos[2];
1065:   if (v) {
1066:     /*
1067:      v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1068:      v[1] = z*cos(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1069:      v[2] = -(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z);
1070:      */
1071:     /*
1072:      v[0] = PetscSinReal(PETSC_PI*x);
1073:      v[1] = PetscSinReal(PETSC_PI*y);
1074:      v[2] = PetscSinReal(PETSC_PI*z);
1075:      */
1076:     v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1077:     v[1] = z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1078:     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;
1079:   }
1080:   if (p) *p = PetscPowRealInt(x,2) + PetscPowRealInt(y,2) + PetscPowRealInt(z,2);
1081:   if (eta) {
1082:     /**eta = PetscExpReal(-theta*(1.0 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));*/
1083:     *eta = 1.0;
1084:   }
1085:   if (Fm) {
1086:     /*
1087:      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))) ;
1088:      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)));
1089:      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)))  ;
1090:      */
1091:     /*
1092:      Fm[0]=-2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*x);
1093:      Fm[1]=-2*y - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*y);
1094:      Fm[2]=-2*z - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*z);
1095:      */
1096:     /*
1097:      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) ;
1098:      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);
1099:      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) ;
1100:      */
1101:     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);
1102:     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);
1103:     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);
1104:   }
1105:   if (Fc) {
1106:     /**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) ;*/
1107:     /**Fc = PETSC_PI*PetscCosReal(PETSC_PI*x) + PETSC_PI*PetscCosReal(PETSC_PI*y) + PETSC_PI*PetscCosReal(PETSC_PI*z);*/
1108:     /**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);*/
1109:     *Fc = 0.0;
1110:   }
1111: }

1113: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx,PetscInt my,PetscInt mz,DM *_da,Vec *_X)
1114: {
1115:   DM             da,cda;
1116:   Vec            X;
1117:   StokesDOF      ***_stokes;
1118:   Vec            coords;
1119:   DMDACoor3d     ***_coords;
1120:   PetscInt       si,sj,sk,ei,ej,ek,i,j,k;

1124:   DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1125:                       mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,4,1,NULL,NULL,NULL,&da);
1126:   DMSetFromOptions(da);
1127:   DMSetUp(da);
1128:   DMDASetFieldName(da,0,"anlytic_Vx");
1129:   DMDASetFieldName(da,1,"anlytic_Vy");
1130:   DMDASetFieldName(da,2,"anlytic_Vz");
1131:   DMDASetFieldName(da,3,"analytic_P");

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

1135:   DMGetCoordinatesLocal(da,&coords);
1136:   DMGetCoordinateDM(da,&cda);
1137:   DMDAVecGetArray(cda,coords,&_coords);

1139:   DMCreateGlobalVector(da,&X);
1140:   DMDAVecGetArray(da,X,&_stokes);

1142:   DMDAGetCorners(da,&si,&sj,&sk,&ei,&ej,&ek);
1143:   for (k = sk; k < sk+ek; k++) {
1144:     for (j = sj; j < sj+ej; j++) {
1145:       for (i = si; i < si+ei; i++) {
1146:         PetscReal pos[NSD],pressure,vel[NSD];

1148:         pos[0] = PetscRealPart(_coords[k][j][i].x);
1149:         pos[1] = PetscRealPart(_coords[k][j][i].y);
1150:         pos[2] = PetscRealPart(_coords[k][j][i].z);

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

1154:         _stokes[k][j][i].u_dof = vel[0];
1155:         _stokes[k][j][i].v_dof = vel[1];
1156:         _stokes[k][j][i].w_dof = vel[2];
1157:         _stokes[k][j][i].p_dof = pressure;
1158:       }
1159:     }
1160:   }
1161:   DMDAVecRestoreArray(da,X,&_stokes);
1162:   DMDAVecRestoreArray(cda,coords,&_coords);

1164:   *_da = da;
1165:   *_X  = X;
1166:   return(0);
1167: }

1169: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da,Vec X,Vec X_analytic)
1170: {
1171:   DM          cda;
1172:   Vec         coords,X_analytic_local,X_local;
1173:   DMDACoor3d  ***_coords;
1174:   PetscInt    sex,sey,sez,mx,my,mz;
1175:   PetscInt    ei,ej,ek;
1176:   PetscScalar el_coords[NODES_PER_EL*NSD];
1177:   StokesDOF   ***stokes_analytic,***stokes;
1178:   StokesDOF   stokes_analytic_e[NODES_PER_EL],stokes_e[NODES_PER_EL];

1180:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1181:   PetscScalar    Ni_p[NODES_PER_EL];
1182:   PetscInt       ngp;
1183:   PetscScalar    gp_xi[GAUSS_POINTS][NSD];
1184:   PetscScalar    gp_weight[GAUSS_POINTS];
1185:   PetscInt       p,i;
1186:   PetscScalar    J_p,fac;
1187:   PetscScalar    h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1188:   PetscScalar    tint_p_ms,tint_p,int_p_ms,int_p;
1189:   PetscInt       M;
1190:   PetscReal      xymin[NSD],xymax[NSD];

1194:   /* define quadrature rule */
1195:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

1197:   /* setup for coords */
1198:   DMGetCoordinateDM(stokes_da,&cda);
1199:   DMGetCoordinatesLocal(stokes_da,&coords);
1200:   DMDAVecGetArray(cda,coords,&_coords);

1202:   /* setup for analytic */
1203:   DMCreateLocalVector(stokes_da,&X_analytic_local);
1204:   DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1205:   DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1206:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1208:   /* setup for solution */
1209:   DMCreateLocalVector(stokes_da,&X_local);
1210:   DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1211:   DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1212:   DMDAVecGetArray(stokes_da,X_local,&stokes);

1214:   DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1215:   DMDAGetBoundingBox(stokes_da,xymin,xymax);

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

1219:   tp_L2     = tu_L2 = tu_H1 = 0.0;
1220:   tint_p_ms = tint_p = 0.0;

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

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:         for (p = 0; p < ngp; p++) {
1234:           ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1235:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1236:           ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1237:           fac = gp_weight[p]*J_p;

1239:           for (i = 0; i < NODES_PER_EL; i++) {
1240:             tint_p_ms = tint_p_ms+fac*Ni_p[i]*stokes_analytic_e[i].p_dof;
1241:             tint_p    = tint_p   +fac*Ni_p[i]*stokes_e[i].p_dof;
1242:           }
1243:         }
1244:       }
1245:     }
1246:   }
1247:   MPI_Allreduce(&tint_p_ms,&int_p_ms,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1248:   MPI_Allreduce(&tint_p,&int_p,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);

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

1252:   /* remove mine and add manufacture one */
1253:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1254:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);

1256:   {
1257:     PetscInt    k,L,dof;
1258:     PetscScalar *fields;
1259:     DMDAGetInfo(stokes_da,0, 0,0,0, 0,0,0, &dof,0,0,0,0,0);

1261:     VecGetLocalSize(X_local,&L);
1262:     VecGetArray(X_local,&fields);
1263:     for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1264:     VecRestoreArray(X_local,&fields);

1266:     VecGetLocalSize(X,&L);
1267:     VecGetArray(X,&fields);
1268:     for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1269:     VecRestoreArray(X,&fields);
1270:   }

1272:   DMDAVecGetArray(stokes_da,X_local,&stokes);
1273:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1275:   for (ek = sez; ek < sez+mz; ek++) {
1276:     for (ej = sey; ej < sey+my; ej++) {
1277:       for (ei = sex; ei < sex+mx; ei++) {
1278:         /* get coords for the element */
1279:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1280:         StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1281:         StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);

1283:         /* evaluate integral */
1284:         p_e_L2 = 0.0;
1285:         u_e_L2 = 0.0;
1286:         u_e_H1 = 0.0;
1287:         for (p = 0; p < ngp; p++) {
1288:           ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1289:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1290:           ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1291:           fac = gp_weight[p]*J_p;

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

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

1298:             u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1299:             v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1300:             w_error = stokes_e[i].w_dof-stokes_analytic_e[i].w_dof;
1301:             u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error+w_error*w_error);

1303:             u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error              /* du/dx */
1304:                                  +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error
1305:                                  +GNx_p[2][i]*u_error*GNx_p[2][i]*u_error
1306:                                  +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error               /* dv/dx */
1307:                                  +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error
1308:                                  +GNx_p[2][i]*v_error*GNx_p[2][i]*v_error
1309:                                  +GNx_p[0][i]*w_error*GNx_p[0][i]*w_error               /* dw/dx */
1310:                                  +GNx_p[1][i]*w_error*GNx_p[1][i]*w_error
1311:                                  +GNx_p[2][i]*w_error*GNx_p[2][i]*w_error);
1312:           }
1313:         }

1315:         tp_L2 += p_e_L2;
1316:         tu_L2 += u_e_L2;
1317:         tu_H1 += u_e_H1;
1318:       }
1319:     }
1320:   }
1321:   MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1322:   MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1323:   MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1324:   p_L2 = PetscSqrtScalar(p_L2);
1325:   u_L2 = PetscSqrtScalar(u_L2);
1326:   u_H1 = PetscSqrtScalar(u_H1);

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


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

1333:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1334:   VecDestroy(&X_analytic_local);
1335:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1336:   VecDestroy(&X_local);
1337:   return(0);
1338: }

1340: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da,Vec FIELD,const char file_prefix[])
1341: {
1342:   char           vtk_filename[PETSC_MAX_PATH_LEN];
1343:   PetscMPIInt    rank;
1344:   MPI_Comm       comm;
1345:   FILE           *vtk_fp = NULL;
1346:   PetscInt       si,sj,sk,nx,ny,nz,i;
1347:   PetscInt       f,n_fields,N;
1348:   DM             cda;
1349:   Vec            coords;
1350:   Vec            l_FIELD;
1351:   PetscScalar    *_L_FIELD;
1352:   PetscInt       memory_offset;
1353:   PetscScalar    *buffer;


1358:   /* create file name */
1359:   PetscObjectGetComm((PetscObject)da,&comm);
1360:   MPI_Comm_rank(comm,&rank);
1361:   PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"subdomain-%s-p%1.4d.vts",file_prefix,rank);

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

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

1369:   /* coords */
1370:   DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1371:   N    = nx * ny * nz;

1373: #if defined(PETSC_WORDS_BIGENDIAN)
1374:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1375: #else
1376:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1377: #endif
1378:   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);
1379:   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);

1381:   memory_offset = 0;

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

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

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

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

1395:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      <PointData Scalars=\" ");
1396:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_fields,0,0,0,0,0);
1397:   for (f=0; f<n_fields; f++) {
1398:     const char *field_name;
1399:     DMDAGetFieldName(da,f,&field_name);
1400:     PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"%s ",field_name);
1401:   }
1402:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\">\n");

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

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

1412:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      </PointData>\n");
1413:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </Piece>\n");
1414:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  </StructuredGrid>\n");

1416:   PetscMalloc1(N,&buffer);
1417:   DMGetLocalVector(da,&l_FIELD);
1418:   DMGlobalToLocalBegin(da, FIELD,INSERT_VALUES,l_FIELD);
1419:   DMGlobalToLocalEnd(da,FIELD,INSERT_VALUES,l_FIELD);
1420:   VecGetArray(l_FIELD,&_L_FIELD);

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

1425:   /* write coordinates */
1426:   {
1427:     int         length = sizeof(PetscScalar)*N*3;
1428:     PetscScalar *allcoords;

1430:     fwrite(&length,sizeof(int),1,vtk_fp);
1431:     VecGetArray(coords,&allcoords);
1432:     fwrite(allcoords,sizeof(PetscScalar),3*N,vtk_fp);
1433:     VecRestoreArray(coords,&allcoords);
1434:   }
1435:   /* write fields */
1436:   for (f=0; f<n_fields; f++) {
1437:     int length = sizeof(PetscScalar)*N;
1438:     fwrite(&length,sizeof(int),1,vtk_fp);
1439:     /* load */
1440:     for (i=0; i<N; i++) buffer[i] = _L_FIELD[n_fields*i + f];

1442:     /* write */
1443:     fwrite(buffer,sizeof(PetscScalar),N,vtk_fp);
1444:   }
1445:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\n  </AppendedData>\n");

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

1449:   PetscFree(buffer);
1450:   VecRestoreArray(l_FIELD,&_L_FIELD);
1451:   DMRestoreLocalVector(da,&l_FIELD);

1453:   if (vtk_fp) {
1454:     fclose(vtk_fp);
1455:     vtk_fp = NULL;
1456:   }

1458:   return(0);
1459: }

1461: PetscErrorCode DAViewVTK_write_PieceExtend(FILE *vtk_fp,PetscInt indent_level,DM da,const char local_file_prefix[])
1462: {
1463:   PetscMPIInt    size,rank;
1464:   MPI_Comm       comm;
1465:   const PetscInt *lx,*ly,*lz;
1466:   PetscInt       M,N,P,pM,pN,pP,sum,*olx,*oly,*olz;
1467:   PetscInt       *osx,*osy,*osz,*oex,*oey,*oez;
1468:   PetscInt       i,j,k,II,stencil;

1472:   /* create file name */
1473:   PetscObjectGetComm((PetscObject)da,&comm);
1474:   MPI_Comm_size(comm,&size);
1475:   MPI_Comm_rank(comm,&rank);

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

1480:   /* generate start,end list */
1481:   PetscMalloc1(pM+1,&olx);
1482:   PetscMalloc1(pN+1,&oly);
1483:   PetscMalloc1(pP+1,&olz);
1484:   sum  = 0;
1485:   for (i=0; i<pM; i++) {
1486:     olx[i] = sum;
1487:     sum    = sum + lx[i];
1488:   }
1489:   olx[pM] = sum;
1490:   sum     = 0;
1491:   for (i=0; i<pN; i++) {
1492:     oly[i] = sum;
1493:     sum    = sum + ly[i];
1494:   }
1495:   oly[pN] = sum;
1496:   sum     = 0;
1497:   for (i=0; i<pP; i++) {
1498:     olz[i] = sum;
1499:     sum    = sum + lz[i];
1500:   }
1501:   olz[pP] = sum;

1503:   PetscMalloc1(pM,&osx);
1504:   PetscMalloc1(pN,&osy);
1505:   PetscMalloc1(pP,&osz);
1506:   PetscMalloc1(pM,&oex);
1507:   PetscMalloc1(pN,&oey);
1508:   PetscMalloc1(pP,&oez);
1509:   for (i=0; i<pM; i++) {
1510:     osx[i] = olx[i] - stencil;
1511:     oex[i] = olx[i] + lx[i] + stencil;
1512:     if (osx[i]<0) osx[i]=0;
1513:     if (oex[i]>M) oex[i]=M;
1514:   }

1516:   for (i=0; i<pN; i++) {
1517:     osy[i] = oly[i] - stencil;
1518:     oey[i] = oly[i] + ly[i] + stencil;
1519:     if (osy[i]<0)osy[i]=0;
1520:     if (oey[i]>M)oey[i]=N;
1521:   }
1522:   for (i=0; i<pP; i++) {
1523:     osz[i] = olz[i] - stencil;
1524:     oez[i] = olz[i] + lz[i] + stencil;
1525:     if (osz[i]<0) osz[i]=0;
1526:     if (oez[i]>P) oez[i]=P;
1527:   }

1529:   for (k=0; k<pP; k++) {
1530:     for (j=0; j<pN; j++) {
1531:       for (i=0; i<pM; i++) {
1532:         char     name[PETSC_MAX_PATH_LEN];
1533:         PetscInt procid = i + j*pM + k*pM*pN; /* convert proc(i,j,k) to pid */
1534:         PetscSNPrintf(name,sizeof(name),"subdomain-%s-p%1.4d.vts",local_file_prefix,procid);
1535:         for (II=0; II<indent_level; II++) PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  ");

1537:         PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<Piece Extent=\"%d %d %d %d %d %d\"      Source=\"%s\"/>\n",
1538:                      osx[i],oex[i]-1,
1539:                      osy[j],oey[j]-1,
1540:                      osz[k],oez[k]-1,name);
1541:       }
1542:     }
1543:   }
1544:   PetscFree(olx);
1545:   PetscFree(oly);
1546:   PetscFree(olz);
1547:   PetscFree(osx);
1548:   PetscFree(osy);
1549:   PetscFree(osz);
1550:   PetscFree(oex);
1551:   PetscFree(oey);
1552:   PetscFree(oez);
1553:   return(0);
1554: }

1556: PetscErrorCode DAView_3DVTK_PStructuredGrid(DM da,const char file_prefix[],const char local_file_prefix[])
1557: {
1558:   MPI_Comm       comm;
1559:   PetscMPIInt    size,rank;
1560:   char           vtk_filename[PETSC_MAX_PATH_LEN];
1561:   FILE           *vtk_fp = NULL;
1562:   PetscInt       M,N,P,si,sj,sk,nx,ny,nz;
1563:   PetscInt       i,dofs;

1567:   /* only master generates this file */
1568:   PetscObjectGetComm((PetscObject)da,&comm);
1569:   MPI_Comm_size(comm,&size);
1570:   MPI_Comm_rank(comm,&rank);

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

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

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

1582: #if defined(PETSC_WORDS_BIGENDIAN)
1583:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1584: #else
1585:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1586: #endif

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

1593:   /* DUMP THE CELL REFERENCES */
1594:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    <PCellData>\n");
1595:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </PCellData>\n");

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

1601:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    <PPointData>\n");
1602:   for (i=0; i<dofs; i++) {
1603:     const char *fieldname;
1604:     DMDAGetFieldName(da,i,&fieldname);
1605:     PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n",fieldname);
1606:   }
1607:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </PPointData>\n");

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

1612:   /* close the file */
1613:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  </PStructuredGrid>\n");
1614:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");

1616:   if (vtk_fp) {
1617:     fclose(vtk_fp);
1618:     vtk_fp = NULL;
1619:   }
1620:   return(0);
1621: }

1623: PetscErrorCode DAView3DPVTS(DM da, Vec x,const char NAME[])
1624: {
1625:   char           vts_filename[PETSC_MAX_PATH_LEN];
1626:   char           pvts_filename[PETSC_MAX_PATH_LEN];

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

1633:   PetscSNPrintf(pvts_filename,sizeof(pvts_filename),"%s-mesh",NAME);
1634:   DAView_3DVTK_PStructuredGrid(da,pvts_filename,vts_filename);
1635:   return(0);
1636: }

1638: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
1639: {
1641:   PetscReal      norms[4];
1642:   Vec            Br,v,w;
1643:   Mat            A;

1646:   KSPGetOperators(ksp,&A,NULL);
1647:   MatCreateVecs(A,&w,&v);

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

1651:   VecStrideNorm(Br,0,NORM_2,&norms[0]);
1652:   VecStrideNorm(Br,1,NORM_2,&norms[1]);
1653:   VecStrideNorm(Br,2,NORM_2,&norms[2]);
1654:   VecStrideNorm(Br,3,NORM_2,&norms[3]);

1656:   VecDestroy(&v);
1657:   VecDestroy(&w);

1659:   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]);
1660:   return(0);
1661: }

1663: static PetscErrorCode PCMGSetupViaCoarsen(PC pc,DM da_fine)
1664: {
1665:   PetscInt              nlevels,k;
1666:   PETSC_UNUSED PetscInt finest;
1667:   DM                    *da_list,*daclist;
1668:   Mat                   R;
1669:   PetscErrorCode        ierr;

1672:   nlevels = 1;
1673:   PetscOptionsGetInt(NULL,NULL,"-levels",&nlevels,0);

1675:   PetscMalloc1(nlevels,&da_list);
1676:   for (k=0; k<nlevels; k++) da_list[k] = NULL;
1677:   PetscMalloc1(nlevels,&daclist);
1678:   for (k=0; k<nlevels; k++) daclist[k] = NULL;

1680:   /* finest grid is nlevels - 1 */
1681:   finest     = nlevels - 1;
1682:   daclist[0] = da_fine;
1683:   PetscObjectReference((PetscObject)da_fine);
1684:   DMCoarsenHierarchy(da_fine,nlevels-1,&daclist[1]);
1685:   for (k=0; k<nlevels; k++) {
1686:     da_list[k] = daclist[nlevels-1-k];
1687:     DMDASetUniformCoordinates(da_list[k],0.0,1.0,0.0,1.0,0.0,1.0);
1688:   }

1690:   PCMGSetLevels(pc,nlevels,NULL);
1691:   PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
1692:   PCMGSetGalerkin(pc,PC_MG_GALERKIN_PMAT);

1694:   for (k=1; k<nlevels; k++) {
1695:     DMCreateInterpolation(da_list[k-1],da_list[k],&R,NULL);
1696:     PCMGSetInterpolation(pc,k,R);
1697:     MatDestroy(&R);
1698:   }

1700:   /* tidy up */
1701:   for (k=0; k<nlevels; k++) {
1702:     DMDestroy(&da_list[k]);
1703:   }
1704:   PetscFree(da_list);
1705:   PetscFree(daclist);
1706:   return(0);
1707: }

1709: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx,PetscInt my,PetscInt mz)
1710: {
1711:   DM             da_Stokes;
1712:   PetscInt       u_dof,p_dof,dof,stencil_width;
1713:   Mat            A,B;
1714:   PetscInt       mxl,myl,mzl;
1715:   DM             vel_cda;
1716:   Vec            vel_coords;
1717:   PetscInt       p;
1718:   Vec            f,X;
1719:   DMDACoor3d     ***_vel_coords;
1720:   PetscInt       its;
1721:   KSP            ksp_S;
1722:   PetscInt       model_definition = 0;
1723:   PetscInt       ei,ej,ek,sex,sey,sez,Imx,Jmy,Kmz;
1724:   CellProperties cell_properties;
1725:   PetscBool      write_output = PETSC_FALSE;

1729:   /* Generate the da for velocity and pressure */
1730:   /* Num nodes in each direction is mx+1, my+1, mz+1 */
1731:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1732:   p_dof         = P_DOFS; /* p - pressure */
1733:   dof           = u_dof+p_dof;
1734:   stencil_width = 1;
1735:   DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1736:                                mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,NULL,&da_Stokes);
1737:   DMSetFromOptions(da_Stokes);
1738:   DMSetUp(da_Stokes);
1739:   DMDASetFieldName(da_Stokes,0,"Vx");
1740:   DMDASetFieldName(da_Stokes,1,"Vy");
1741:   DMDASetFieldName(da_Stokes,2,"Vz");
1742:   DMDASetFieldName(da_Stokes,3,"P");

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

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

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

1753:   /* interpolate the coordinates to quadrature points */
1754:   DMGetCoordinateDM(da_Stokes,&vel_cda);
1755:   DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1756:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1757:   DMDAGetElementCorners(da_Stokes,&sex,&sey,&sez,&Imx,&Jmy,&Kmz);
1758:   for (ek = sez; ek < sez+Kmz; ek++) {
1759:     for (ej = sey; ej < sey+Jmy; ej++) {
1760:       for (ei = sex; ei < sex+Imx; ei++) {
1761:         /* get coords for the element */
1762:         PetscInt               ngp;
1763:         PetscScalar            gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
1764:         PetscScalar            el_coords[NSD*NODES_PER_EL];
1765:         GaussPointCoefficients *cell;

1767:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1768:         GetElementCoords3D(_vel_coords,ei,ej,ek,el_coords);
1769:         ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

1771:         for (p = 0; p < GAUSS_POINTS; p++) {
1772:           PetscScalar xi_p[NSD],Ni_p[NODES_PER_EL];
1773:           PetscScalar gp_x,gp_y,gp_z;
1774:           PetscInt    n;

1776:           xi_p[0] = gp_xi[p][0];
1777:           xi_p[1] = gp_xi[p][1];
1778:           xi_p[2] = gp_xi[p][2];
1779:           ShapeFunctionQ13D_Evaluate(xi_p,Ni_p);

1781:           gp_x = gp_y = gp_z = 0.0;
1782:           for (n = 0; n < NODES_PER_EL; n++) {
1783:             gp_x = gp_x+Ni_p[n]*el_coords[NSD*n];
1784:             gp_y = gp_y+Ni_p[n]*el_coords[NSD*n+1];
1785:             gp_z = gp_z+Ni_p[n]*el_coords[NSD*n+2];
1786:           }
1787:           cell->gp_coords[NSD*p]   = gp_x;
1788:           cell->gp_coords[NSD*p+1] = gp_y;
1789:           cell->gp_coords[NSD*p+2] = gp_z;
1790:         }
1791:       }
1792:     }
1793:   }

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

1797:   switch (model_definition) {
1798:   case 0: /* isoviscous */
1799:     for (ek = sez; ek < sez+Kmz; ek++) {
1800:       for (ej = sey; ej < sey+Jmy; ej++) {
1801:         for (ei = sex; ei < sex+Imx; ei++) {
1802:           GaussPointCoefficients *cell;

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

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

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

1822:   case 1: /* manufactured */
1823:     for (ek = sez; ek < sez+Kmz; ek++) {
1824:       for (ej = sey; ej < sey+Jmy; ej++) {
1825:         for (ei = sex; ei < sex+Imx; ei++) {
1826:           PetscReal              eta,Fm[NSD],Fc,pos[NSD];
1827:           GaussPointCoefficients *cell;

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

1835:             pos[0] = coord_x;
1836:             pos[1] = coord_y;
1837:             pos[2] = coord_z;

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

1842:             cell->fx[p] = Fm[0];
1843:             cell->fy[p] = Fm[1];
1844:             cell->fz[p] = Fm[2];
1845:             cell->hc[p] = Fc;
1846:           }
1847:         }
1848:       }
1849:     }
1850:     break;

1852:   case 2: /* solcx */
1853:     for (ek = sez; ek < sez+Kmz; ek++) {
1854:       for (ej = sey; ej < sey+Jmy; ej++) {
1855:         for (ei = sex; ei < sex+Imx; ei++) {
1856:           GaussPointCoefficients *cell;

1858:           CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1859:           for (p = 0; p < GAUSS_POINTS; p++) {
1860:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1861:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1862:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);

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

1866:             cell->fx[p] = 0.0;
1867:             cell->fy[p] = -PetscSinReal(3.0*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1868:             cell->fz[p] = 0.0*coord_z;
1869:             cell->hc[p] = 0.0;
1870:           }
1871:         }
1872:       }
1873:     }
1874:     break;

1876:   case 3: /* sinker */
1877:     for (ek = sez; ek < sez+Kmz; ek++) {
1878:       for (ej = sey; ej < sey+Jmy; ej++) {
1879:         for (ei = sex; ei < sex+Imx; ei++) {
1880:           GaussPointCoefficients *cell;

1882:           CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1883:           for (p = 0; p < GAUSS_POINTS; p++) {
1884:             PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p]);
1885:             PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);
1886:             PetscReal zp = PetscRealPart(cell->gp_coords[NSD*p+2]);

1888:             cell->eta[p] = 1.0e-2;
1889:             cell->fx[p]  = 0.0;
1890:             cell->fy[p]  = 0.0;
1891:             cell->fz[p]  = 0.0;
1892:             cell->hc[p]  = 0.0;

1894:             if ((PetscAbsReal(xp-0.5) < 0.2) &&
1895:                 (PetscAbsReal(yp-0.5) < 0.2) &&
1896:                 (PetscAbsReal(zp-0.5) < 0.2)) {
1897:               cell->eta[p] = 1.0;
1898:               cell->fz[p]  = 1.0;
1899:             }

1901:           }
1902:         }
1903:       }
1904:     }
1905:     break;

1907:   default:
1908:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"No default model is supported. Choose either -model {0,1,2,3}");
1909:     break;
1910:   }

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

1914:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1915:   DMSetMatType(da_Stokes,MATAIJ);
1916:   DMCreateMatrix(da_Stokes,&A);
1917:   DMCreateMatrix(da_Stokes,&B);
1918:   DMCreateGlobalVector(da_Stokes,&X);
1919:   DMCreateGlobalVector(da_Stokes,&f);

1921:   /* assemble A11 */
1922:   MatZeroEntries(A);
1923:   MatZeroEntries(B);
1924:   VecZeroEntries(f);

1926:   AssembleA_Stokes(A,da_Stokes,cell_properties);
1927:   AssembleA_PCStokes(B,da_Stokes,cell_properties);
1928:   /* build force vector */
1929:   AssembleF_Stokes(f,da_Stokes,cell_properties);

1931:   /* SOLVE */
1932:   KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1933:   KSPSetOptionsPrefix(ksp_S,"stokes_"); /* stokes */
1934:   KSPSetOperators(ksp_S,A,B);
1935:   KSPSetFromOptions(ksp_S);

1937:   {
1938:     PC             pc;
1939:     const PetscInt ufields[] = {0,1,2},pfields[] = {3};
1940:     KSPGetPC(ksp_S,&pc);
1941:     PCFieldSplitSetBlockSize(pc,4);
1942:     PCFieldSplitSetFields(pc,"u",3,ufields,ufields);
1943:     PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1944:   }

1946:   {
1947:     PC        pc;
1948:     PetscBool same = PETSC_FALSE;
1949:     KSPGetPC(ksp_S,&pc);
1950:     PetscObjectTypeCompare((PetscObject)pc,PCMG,&same);
1951:     if (same) {
1952:       PCMGSetupViaCoarsen(pc,da_Stokes);
1953:     }
1954:   }

1956:   {
1957:     PetscBool stokes_monitor = PETSC_FALSE;
1958:     PetscOptionsGetBool(NULL,NULL,"-stokes_ksp_monitor_blocks",&stokes_monitor,0);
1959:     if (stokes_monitor) {
1960:       KSPMonitorSet(ksp_S,KSPMonitorStokesBlocks,NULL,NULL);
1961:     }
1962:   }
1963:   KSPSolve(ksp_S,f,X);

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

1981:   /* verify */
1982:   if (model_definition == 1) {
1983:     DM  da_Stokes_analytic;
1984:     Vec X_analytic;

1986:     DMDACreateManufacturedSolution(mx,my,mz,&da_Stokes_analytic,&X_analytic);
1987:     if (write_output) {
1988:       DAView3DPVTS(da_Stokes_analytic,X_analytic,"ms");
1989:     }
1990:     DMDAIntegrateErrors3D(da_Stokes_analytic,X,X_analytic);
1991:     if (write_output) {
1992:       DAView3DPVTS(da_Stokes,X,"up2");
1993:     }
1994:     DMDestroy(&da_Stokes_analytic);
1995:     VecDestroy(&X_analytic);
1996:   }

1998:   KSPDestroy(&ksp_S);
1999:   VecDestroy(&X);
2000:   VecDestroy(&f);
2001:   MatDestroy(&A);
2002:   MatDestroy(&B);

2004:   CellPropertiesDestroy(&cell_properties);
2005:   DMDestroy(&da_Stokes);
2006:   return(0);
2007: }

2009: int main(int argc,char **args)
2010: {
2012:   PetscInt       mx,my,mz;

2014:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
2015:   mx   = my = mz = 10;
2016:   PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
2017:   my   = mx; mz = mx;
2018:   PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
2019:   PetscOptionsGetInt(NULL,NULL,"-mz",&mz,NULL);
2020:   solve_stokes_3d_coupled(mx,my,mz);
2021:   PetscFinalize();
2022:   return ierr;
2023: }