Actual source code: ex42.c

petsc-3.4.5 2014-06-29
  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
 14: #include <petsctime.h>

 16: #define PROFILE_TIMING
 17: #define ASSEMBLE_LOWER_TRIANGULAR

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

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

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

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

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

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

 62:   PetscMalloc(sizeof(struct _p_CellProperties),&cells);

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

 66:   cells->mx     = mx;
 67:   cells->my     = my;
 68:   cells->mz     = mz;
 69:   cells->ncells = mx * my * mz;
 70:   cells->sex    = sex;
 71:   cells->sey    = sey;
 72:   cells->sez    = sez;

 74:   PetscMalloc(sizeof(GaussPointCoefficients)*mx*my*mz,&cells->gpc);

 76:   *C = cells;
 77:   return(0);
 78: }

 82: PetscErrorCode CellPropertiesDestroy(CellProperties *C)
 83: {
 85:   CellProperties cells;

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

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

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

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

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

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

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

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

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

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

171:   t4  = A[2][0] * A[0][1];
172:   t6  = A[2][0] * A[0][2];
173:   t8  = A[1][0] * A[0][1];
174:   t10 = A[1][0] * A[0][2];
175:   t12 = A[0][0] * A[1][1];
176:   t14 = A[0][0] * A[1][2];
177:   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]);

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

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

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

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

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

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

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

222:   matrix_inverse_3x3(JJ,iJ);

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

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

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

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

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

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

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

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

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

287: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
288: {
289:   PetscInt       si,sj,sk;

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

295:   if (sx) {
296:     *sx = si;
297:     if (si != 0) *sx = si+1;
298:   }
299:   if (sy) {
300:     *sy = sj;
301:     if (sj != 0) *sy = sj+1;
302:   }
303:   if (sz) {
304:     *sz = sk;
305:     if (sk != 0) *sz = sk+1;
306:   }
307:   DMDAGetLocalElementSize(da,mx,my,mz);
308:   return(0);
309: }

311: /*
312:  i,j are the element indices
313:  The unknown is a vector quantity.
314:  The s[].c is used to indicate the degree of freedom.
315:  */
318: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j,PetscInt k)
319: {
320:   PetscInt n;

323:   /* velocity */
324:   n = 0;
325:   /* node 0 */
326:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++; /* Vx0 */
327:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++; /* Vy0 */
328:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++; /* Vz0 */

330:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
331:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
332:   s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;

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

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

342:   /* */
343:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++; /* Vx4 */
344:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++; /* Vy4 */
345:   s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++; /* Vz4 */

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

351:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
352:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
353:   s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;

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

359:   /* pressure */
360:   n = 0;

362:   s_p[n].i = i;   s_p[n].j = j;   s_p[n].k = k; s_p[n].c = 3; n++; /* P0 */
363:   s_p[n].i = i;   s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
364:   s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
365:   s_p[n].i = i+1; s_p[n].j = j;   s_p[n].k = k; s_p[n].c = 3; n++;

367:   s_p[n].i = i;   s_p[n].j = j;   s_p[n].k = k+1; s_p[n].c = 3; n++; /* P0 */
368:   s_p[n].i = i;   s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
369:   s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
370:   s_p[n].i = i+1; s_p[n].j = j;   s_p[n].k = k+1; s_p[n].c = 3; n++;
371:   return(0);
372: }

376: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords,PetscInt i,PetscInt j,PetscInt k,PetscScalar el_coord[])
377: {
379:   /* get coords for the element */
380:   el_coord[0] = coords[k][j][i].x;
381:   el_coord[1] = coords[k][j][i].y;
382:   el_coord[2] = coords[k][j][i].z;

384:   el_coord[3] = coords[k][j+1][i].x;
385:   el_coord[4] = coords[k][j+1][i].y;
386:   el_coord[5] = coords[k][j+1][i].z;

388:   el_coord[6] = coords[k][j+1][i+1].x;
389:   el_coord[7] = coords[k][j+1][i+1].y;
390:   el_coord[8] = coords[k][j+1][i+1].z;

392:   el_coord[9]  = coords[k][j][i+1].x;
393:   el_coord[10] = coords[k][j][i+1].y;
394:   el_coord[11] = coords[k][j][i+1].z;

396:   el_coord[12] = coords[k+1][j][i].x;
397:   el_coord[13] = coords[k+1][j][i].y;
398:   el_coord[14] = coords[k+1][j][i].z;

400:   el_coord[15] = coords[k+1][j+1][i].x;
401:   el_coord[16] = coords[k+1][j+1][i].y;
402:   el_coord[17] = coords[k+1][j+1][i].z;

404:   el_coord[18] = coords[k+1][j+1][i+1].x;
405:   el_coord[19] = coords[k+1][j+1][i+1].y;
406:   el_coord[20] = coords[k+1][j+1][i+1].z;

408:   el_coord[21] = coords[k+1][j][i+1].x;
409:   el_coord[22] = coords[k+1][j][i+1].y;
410:   el_coord[23] = coords[k+1][j][i+1].z;
411:   return(0);
412: }

416: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field,PetscInt i,PetscInt j,PetscInt k,StokesDOF nodal_fields[])
417: {
419:   /* get the nodal fields for u */
420:   nodal_fields[0].u_dof = field[k][j][i].u_dof;
421:   nodal_fields[0].v_dof = field[k][j][i].v_dof;
422:   nodal_fields[0].w_dof = field[k][j][i].w_dof;

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

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

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

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

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

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

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

452:   /* get the nodal fields for p */
453:   nodal_fields[0].p_dof = field[k][j][i].p_dof;
454:   nodal_fields[1].p_dof = field[k][j+1][i].p_dof;
455:   nodal_fields[2].p_dof = field[k][j+1][i+1].p_dof;
456:   nodal_fields[3].p_dof = field[k][j][i+1].p_dof;

458:   nodal_fields[4].p_dof = field[k+1][j][i].p_dof;
459:   nodal_fields[5].p_dof = field[k+1][j+1][i].p_dof;
460:   nodal_fields[6].p_dof = field[k+1][j+1][i+1].p_dof;
461:   nodal_fields[7].p_dof = field[k+1][j][i+1].p_dof;
462:   return(0);
463: }

465: 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)
466: {
467:   PetscInt              ij;
468:   PETSC_UNUSED PetscInt r,c,nr,nc;

470:   nr = w_NPE*w_dof;
471:   nc = u_NPE*u_dof;

473:   r = w_dof*wi+wd;
474:   c = u_dof*ui+ud;

476:   ij = r*nc+c;

478:   return ij;
479: }

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

488:   for (n = 0; n<NODES_PER_EL; n++) {
489:     I = u_eqn[NSD*n].i;
490:     J = u_eqn[NSD*n].j;
491:     K = u_eqn[NSD*n].k;

493:     fields_F[K][J][I].u_dof = fields_F[K][J][I].u_dof+Fe_u[NSD*n];

495:     I = u_eqn[NSD*n+1].i;
496:     J = u_eqn[NSD*n+1].j;
497:     K = u_eqn[NSD*n+1].k;

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

501:     I = u_eqn[NSD*n+2].i;
502:     J = u_eqn[NSD*n+2].j;
503:     K = u_eqn[NSD*n+2].k;
504:     fields_F[K][J][I].w_dof = fields_F[K][J][I].w_dof+Fe_u[NSD*n+2];

506:     I = p_eqn[n].i;
507:     J = p_eqn[n].j;
508:     K = p_eqn[n].k;

510:     fields_F[K][J][I].p_dof = fields_F[K][J][I].p_dof+Fe_p[n];

512:   }
513:   return(0);
514: }

516: static void FormStressOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
517: {
518:   PetscInt       ngp;
519:   PetscScalar    gp_xi[GAUSS_POINTS][NSD];
520:   PetscScalar    gp_weight[GAUSS_POINTS];
521:   PetscInt       p,i,j,k;
522:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
523:   PetscScalar    J_p,tildeD[6];
524:   PetscScalar    B[6][U_DOFS*NODES_PER_EL];
525:   const PetscInt nvdof = U_DOFS*NODES_PER_EL;

527:   /* define quadrature rule */
528:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

530:   /* evaluate integral */
531:   for (p = 0; p < ngp; p++) {
532:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
533:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);

535:     for (i = 0; i < NODES_PER_EL; i++) {
536:       PetscScalar d_dx_i = GNx_p[0][i];
537:       PetscScalar d_dy_i = GNx_p[1][i];
538:       PetscScalar d_dz_i = GNx_p[2][i];

540:       B[0][3*i] = d_dx_i; B[0][3*i+1] = 0.0;     B[0][3*i+2] = 0.0;
541:       B[1][3*i] = 0.0;    B[1][3*i+1] = d_dy_i;  B[1][3*i+2] = 0.0;
542:       B[2][3*i] = 0.0;    B[2][3*i+1] = 0.0;     B[2][3*i+2] = d_dz_i;

544:       B[3][3*i] = d_dy_i; B[3][3*i+1] = d_dx_i;  B[3][3*i+2] = 0.0;   /* e_xy */
545:       B[4][3*i] = d_dz_i; B[4][3*i+1] = 0.0;     B[4][3*i+2] = d_dx_i; /* e_xz */
546:       B[5][3*i] = 0.0;    B[5][3*i+1] = d_dz_i;  B[5][3*i+2] = d_dy_i; /* e_yz */
547:     }


550:     tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
551:     tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
552:     tildeD[2] = 2.0*gp_weight[p]*J_p*eta[p];

554:     tildeD[3] =     gp_weight[p]*J_p*eta[p];
555:     tildeD[4] =     gp_weight[p]*J_p*eta[p];
556:     tildeD[5] =     gp_weight[p]*J_p*eta[p];

558:     /* form Bt tildeD B */
559:     /*
560:      Ke_ij = Bt_ik . D_kl . B_lj
561:      = B_ki . D_kl . B_lj
562:      = B_ki . D_kk . B_kj
563:      */
564:     for (i = 0; i < nvdof; i++) {
565:       for (j = i; j < nvdof; j++) {
566:         for (k = 0; k < 6; k++) {
567:           Ke[i*nvdof+j] += B[k][i]*tildeD[k]*B[k][j];
568:         }
569:       }
570:     }

572:   }
573:   /* fill lower triangular part */
574: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
575:   for (i = 0; i < nvdof; i++) {
576:     for (j = i; j < nvdof; j++) {
577:       Ke[j*nvdof+i] = Ke[i*nvdof+j];
578:     }
579:   }
580: #endif
581: }

583: static void FormGradientOperatorQ13D(PetscScalar Ke[],PetscScalar coords[])
584: {
585:   PetscInt    ngp;
586:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
587:   PetscScalar gp_weight[GAUSS_POINTS];
588:   PetscInt    p,i,j,di;
589:   PetscScalar Ni_p[NODES_PER_EL];
590:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
591:   PetscScalar J_p,fac;

593:   /* define quadrature rule */
594:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

596:   /* evaluate integral */
597:   for (p = 0; p < ngp; p++) {
598:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
599:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
600:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
601:     fac = gp_weight[p]*J_p;

603:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
604:       for (di = 0; di < NSD; di++) { /* u dofs */
605:         for (j = 0; j < NODES_PER_EL; j++) {  /* p nodes, p dofs = 1 (ie no loop) */
606:           PetscInt IJ;
607:           IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,3,j,0,NODES_PER_EL,1);

609:           Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
610:         }
611:       }
612:     }
613:   }
614: }

616: static void FormDivergenceOperatorQ13D(PetscScalar De[],PetscScalar coords[])
617: {
618:   PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
619:   PetscInt    i,j;
620:   PetscInt    nr_g,nc_g;

622:   PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
623:   FormGradientOperatorQ13D(Ge,coords);

625:   nr_g = U_DOFS*NODES_PER_EL;
626:   nc_g = P_DOFS*NODES_PER_EL;

628:   for (i = 0; i < nr_g; i++) {
629:     for (j = 0; j < nc_g; j++) {
630:       De[nr_g*j+i] = Ge[nc_g*i+j];
631:     }
632:   }
633: }

635: static void FormStabilisationOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
636: {
637:   PetscInt    ngp;
638:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
639:   PetscScalar gp_weight[GAUSS_POINTS];
640:   PetscInt    p,i,j;
641:   PetscScalar Ni_p[NODES_PER_EL];
642:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
643:   PetscScalar J_p,fac,eta_avg;

645:   /* define quadrature rule */
646:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

648:   /* evaluate integral */
649:   for (p = 0; p < ngp; p++) {
650:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
651:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
652:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
653:     fac = gp_weight[p]*J_p;
654:     /*
655:      for (i = 0; i < NODES_PER_EL; i++) {
656:      for (j = i; j < NODES_PER_EL; j++) {
657:      Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
658:      }
659:      }
660:      */

662:     for (i = 0; i < NODES_PER_EL; i++) {
663:       for (j = 0; j < NODES_PER_EL; j++) {
664:         Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
665:       }
666:     }
667:   }

669:   /* scale */
670:   eta_avg = 0.0;
671:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
672:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
673:   fac     = 1.0/eta_avg;

675:   /*
676:    for (i = 0; i < NODES_PER_EL; i++) {
677:    for (j = i; j < NODES_PER_EL; j++) {
678:    Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
679:    #if defined(ASSEMBLE_LOWER_TRIANGULAR)
680:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
681:    #endif
682:    }
683:    }
684:    */

686:   for (i = 0; i < NODES_PER_EL; i++) {
687:     for (j = 0; j < NODES_PER_EL; j++) {
688:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
689:     }
690:   }
691: }

693: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
694: {
695:   PetscInt    ngp;
696:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
697:   PetscScalar gp_weight[GAUSS_POINTS];
698:   PetscInt    p,i,j;
699:   PetscScalar Ni_p[NODES_PER_EL];
700:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
701:   PetscScalar J_p,fac,eta_avg;

703:   /* define quadrature rule */
704:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

706:   /* evaluate integral */
707:   for (p = 0; p < ngp; p++) {
708:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
709:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
710:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
711:     fac = gp_weight[p]*J_p;

713:     /*
714:      for (i = 0; i < NODES_PER_EL; i++) {
715:      for (j = i; j < NODES_PER_EL; j++) {
716:      Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
717:      }
718:      }
719:      */

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

728:   /* scale */
729:   eta_avg = 0.0;
730:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
731:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
732:   fac     = 1.0/eta_avg;
733:   /*
734:    for (i = 0; i < NODES_PER_EL; i++) {
735:    for (j = i; j < NODES_PER_EL; j++) {
736:    Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
737:    #if defined(ASSEMBLE_LOWER_TRIANGULAR)
738:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
739:    #endif
740:    }
741:    }
742:    */

744:   for (i = 0; i < NODES_PER_EL; i++) {
745:     for (j = 0; j < NODES_PER_EL; j++) {
746:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
747:     }
748:   }
749: }

751: static void FormMomentumRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[],PetscScalar fz[])
752: {
753:   PetscInt    ngp;
754:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
755:   PetscScalar gp_weight[GAUSS_POINTS];
756:   PetscInt    p,i;
757:   PetscScalar Ni_p[NODES_PER_EL];
758:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
759:   PetscScalar J_p,fac;

761:   /* define quadrature rule */
762:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

764:   /* evaluate integral */
765:   for (p = 0; p < ngp; p++) {
766:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
767:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
768:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
769:     fac = gp_weight[p]*J_p;

771:     for (i = 0; i < NODES_PER_EL; i++) {
772:       Fe[NSD*i]   -= fac*Ni_p[i]*fx[p];
773:       Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
774:       Fe[NSD*i+2] -= fac*Ni_p[i]*fz[p];
775:     }
776:   }
777: }

779: static void FormContinuityRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar hc[])
780: {
781:   PetscInt    ngp;
782:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
783:   PetscScalar gp_weight[GAUSS_POINTS];
784:   PetscInt    p,i;
785:   PetscScalar Ni_p[NODES_PER_EL];
786:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
787:   PetscScalar J_p,fac;

789:   /* define quadrature rule */
790:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

792:   /* evaluate integral */
793:   for (p = 0; p < ngp; p++) {
794:     ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
795:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
796:     ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
797:     fac = gp_weight[p]*J_p;

799:     for (i = 0; i < NODES_PER_EL; i++) Fe[i] -= fac*Ni_p[i]*hc[p];
800:   }
801: }

803: #define _ZERO_ROWCOL_i(A,i) {                   \
804:     PetscInt    KK;                             \
805:     PetscScalar tmp = A[24*(i)+(i)];            \
806:     for (KK=0;KK<24;KK++) A[24*(i)+KK]=0.0;     \
807:     for (KK=0;KK<24;KK++) A[24*KK+(i)]=0.0;     \
808:     A[24*(i)+(i)] = tmp;}                       \

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

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

820: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,CellProperties cell_properties)
821: {
822:   DM                     cda;
823:   Vec                    coords;
824:   DMDACoor3d             ***_coords;
825:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS];
826:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS];
827:   PetscInt               sex,sey,sez,mx,my,mz;
828:   PetscInt               ei,ej,ek;
829:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
830:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
831:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
832:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
833:   PetscScalar            el_coords[NODES_PER_EL*NSD];
834:   GaussPointCoefficients *props;
835:   PetscScalar            *prop_eta;
836:   PetscInt               n,M,N,P;
837:   PetscLogDouble         t0,t1;
838:   PetscErrorCode         ierr;

841:   DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
842:   /* setup for coords */
843:   DMGetCoordinateDM(stokes_da,&cda);
844:   DMGetCoordinatesLocal(stokes_da,&coords);
845:   DMDAVecGetArray(cda,coords,&_coords);

847:   DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
848:   PetscTime(&t0);
849:   for (ek = sez; ek < sez+mz; ek++) {
850:     for (ej = sey; ej < sey+my; ej++) {
851:       for (ei = sex; ei < sex+mx; ei++) {
852:         /* get coords for the element */
853:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
854:         /* get cell properties */
855:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
856:         /* get coefficients for the element */
857:         prop_eta = props->eta;

859:         /* initialise element stiffness matrix */
860:         PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
861:         PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
862:         PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
863:         PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);

865:         /* form element stiffness matrix */
866:         FormStressOperatorQ13D(Ae,el_coords,prop_eta);
867:         FormGradientOperatorQ13D(Ge,el_coords);
868:         /*#if defined(ASSEMBLE_LOWER_TRIANGULAR)*/
869:         FormDivergenceOperatorQ13D(De,el_coords);
870:         /*#endif*/
871:         FormStabilisationOperatorQ13D(Ce,el_coords,prop_eta);

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

876:         for (n=0; n<NODES_PER_EL; n++) {
877:           if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
878:             _ZERO_ROWCOL_i(Ae,3*n);
879:             _ZERO_ROW_i   (Ge,3*n);
880:             _ZERO_COL_i   (De,3*n);
881:           }

883:           if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
884:             _ZERO_ROWCOL_i(Ae,3*n+1);
885:             _ZERO_ROW_i   (Ge,3*n+1);
886:             _ZERO_COL_i   (De,3*n+1);
887:           }

889:           if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
890:             _ZERO_ROWCOL_i(Ae,3*n+2);
891:             _ZERO_ROW_i   (Ge,3*n+2);
892:             _ZERO_COL_i   (De,3*n+2);
893:           }
894:         }
895:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
896:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
897:         MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
898:         MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
899:       }
900:     }
901:   }
902:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
903:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

907:   PetscTime(&t1);
908:   return(0);
909: }

913: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,CellProperties cell_properties)
914: {
915:   DM                     cda;
916:   Vec                    coords;
917:   DMDACoor3d             ***_coords;
918:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS];
919:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS];
920:   PetscInt               sex,sey,sez,mx,my,mz;
921:   PetscInt               ei,ej,ek;
922:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
923:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
924:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
925:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
926:   PetscScalar            el_coords[NODES_PER_EL*NSD];
927:   GaussPointCoefficients *props;
928:   PetscScalar            *prop_eta;
929:   PetscInt               n,M,N,P;
930:   PetscErrorCode         ierr;

933:   DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
934:   /* setup for coords */
935:   DMGetCoordinateDM(stokes_da,&cda);
936:   DMGetCoordinatesLocal(stokes_da,&coords);
937:   DMDAVecGetArray(cda,coords,&_coords);

939:   DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
940:   for (ek = sez; ek < sez+mz; ek++) {
941:     for (ej = sey; ej < sey+my; ej++) {
942:       for (ei = sex; ei < sex+mx; ei++) {
943:         /* get coords for the element */
944:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
945:         /* get cell properties */
946:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
947:         /* get coefficients for the element */
948:         prop_eta = props->eta;

950:         /* initialise element stiffness matrix */
951:         PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
952:         PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
953:         PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
954:         PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);

956:         /* form element stiffness matrix */
957:         FormStressOperatorQ13D(Ae,el_coords,prop_eta);
958:         FormGradientOperatorQ13D(Ge,el_coords);
959:         /* FormDivergenceOperatorQ13D(De,el_coords); */
960:         FormScaledMassMatrixOperatorQ13D(Ce,el_coords,prop_eta);

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

965:         for (n=0; n<NODES_PER_EL; n++) {
966:           if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
967:             _ZERO_ROWCOL_i(Ae,3*n);
968:             _ZERO_ROW_i   (Ge,3*n);
969:             _ZERO_COL_i   (De,3*n);
970:           }

972:           if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
973:             _ZERO_ROWCOL_i(Ae,3*n+1);
974:             _ZERO_ROW_i   (Ge,3*n+1);
975:             _ZERO_COL_i   (De,3*n+1);
976:           }

978:           if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
979:             _ZERO_ROWCOL_i(Ae,3*n+2);
980:             _ZERO_ROW_i   (Ge,3*n+2);
981:             _ZERO_COL_i   (De,3*n+2);
982:           }
983:         }
984:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
985:         MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
986:         /*MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);*/
987:         MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
988:       }
989:     }
990:   }
991:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
992:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

1000: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,CellProperties cell_properties)
1001: {
1002:   DM                     cda;
1003:   Vec                    coords;
1004:   DMDACoor3d             ***_coords;
1005:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS];
1006:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS];
1007:   PetscInt               sex,sey,sez,mx,my,mz;
1008:   PetscInt               ei,ej,ek;
1009:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
1010:   PetscScalar            He[NODES_PER_EL*P_DOFS];
1011:   PetscScalar            el_coords[NODES_PER_EL*NSD];
1012:   GaussPointCoefficients *props;
1013:   PetscScalar            *prop_fx,*prop_fy,*prop_fz,*prop_hc;
1014:   Vec                    local_F;
1015:   StokesDOF              ***ff;
1016:   PetscInt               n,M,N,P;
1017:   PetscErrorCode         ierr;

1020:   DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
1021:   /* setup for coords */
1022:   DMGetCoordinateDM(stokes_da,&cda);
1023:   DMGetCoordinatesLocal(stokes_da,&coords);
1024:   DMDAVecGetArray(cda,coords,&_coords);

1026:   /* get acces to the vector */
1027:   DMGetLocalVector(stokes_da,&local_F);
1028:   VecZeroEntries(local_F);
1029:   DMDAVecGetArray(stokes_da,local_F,&ff);

1031:   DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
1032:   for (ek = sez; ek < sez+mz; ek++) {
1033:     for (ej = sey; ej < sey+my; ej++) {
1034:       for (ei = sex; ei < sex+mx; ei++) {
1035:         /* get coords for the element */
1036:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1037:         /* get cell properties */
1038:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
1039:         /* get coefficients for the element */
1040:         prop_fx = props->fx;
1041:         prop_fy = props->fy;
1042:         prop_fz = props->fz;
1043:         prop_hc = props->hc;

1045:         /* initialise element stiffness matrix */
1046:         PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
1047:         PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);

1049:         /* form element stiffness matrix */
1050:         FormMomentumRhsQ13D(Fe,el_coords,prop_fx,prop_fy,prop_fz);
1051:         FormContinuityRhsQ13D(He,el_coords,prop_hc);

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

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

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

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

1064:         DMDASetValuesLocalStencil3D_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
1065:       }
1066:     }
1067:   }
1068:   DMDAVecRestoreArray(stokes_da,local_F,&ff);
1069:   DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
1070:   DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
1071:   DMRestoreLocalVector(stokes_da,&local_F);

1073:   DMDAVecRestoreArray(cda,coords,&_coords);
1074:   return(0);
1075: }

1077: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta,PetscReal *MX,PetscReal *MY,PetscReal *MZ)
1078: {
1079:   *theta = 0.0;
1080:   *MX    = 2.0 * M_PI;
1081:   *MY    = 2.0 * M_PI;
1082:   *MZ    = 2.0 * M_PI;
1083: }
1084: static void evaluate_MS_FrankKamentski(PetscReal pos[],PetscReal v[],PetscReal *p,PetscReal *eta,PetscReal Fm[],PetscReal *Fc)
1085: {
1086:   PetscReal x,y,z;
1087:   PetscReal theta,MX,MY,MZ;

1089:   evaluate_MS_FrankKamentski_constants(&theta,&MX,&MY,&MZ);
1090:   x = pos[0];
1091:   y = pos[1];
1092:   z = pos[2];
1093:   if (v) {
1094:     /*
1095:      v[0] = pow(z,3)*exp(y)*sin(M_PI*x);
1096:      v[1] = z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y);
1097:      v[2] = -(pow(x,3) + pow(y,3))*sin(2.0*M_PI*z);
1098:      */
1099:     /*
1100:      v[0] = sin(M_PI*x);
1101:      v[1] = sin(M_PI*y);
1102:      v[2] = sin(M_PI*z);
1103:      */
1104:     v[0] = PetscPowRealInt(z,3)*exp(y)*sin(M_PI*x);
1105:     v[1] = z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y);
1106:     v[2] = PetscPowRealInt(z,2)*(cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 - M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2) - M_PI*PetscPowRealInt(z,4)*cos(M_PI*x)*exp(y)/4;
1107:   }
1108:   if (p) *p = PetscPowRealInt(x,2) + PetscPowRealInt(y,2) + PetscPowRealInt(z,2);
1109:   if (eta) {
1110:     /**eta = exp(-theta*(1.0 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)));*/
1111:     *eta = 1.0;
1112:   }
1113:   if (Fm) {
1114:     /*
1115:      Fm[0] = -2*x - 2.0*pow(M_PI,2)*pow(z,3)*exp(y)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(M_PI*x) - 0.2*MZ*theta*(-1.5*pow(x,2)*sin(2.0*M_PI*z) + 1.5*pow(z,2)*exp(y)*sin(M_PI*x))*cos(MX*x)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MY*y)*sin(MZ*z) - 0.2*M_PI*MX*theta*pow(z,3)*cos(M_PI*x)*cos(MZ*z)*exp(y)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MX*x)*sin(MY*y) + 2.0*(3.0*z*exp(y)*sin(M_PI*x) - 3.0*M_PI*pow(x,2)*cos(2.0*M_PI*z))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*(0.5*pow(z,3)*exp(y)*sin(M_PI*x) + M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x) - 1.0*z*pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*theta*(1 + 0.1*MY*cos(MX*x)*cos(MY*y)*cos(MZ*z))*(0.5*pow(z,3)*exp(y)*sin(M_PI*x) - 1.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) ;
1116:      Fm[1] = -2*y - 0.2*MX*theta*(0.5*pow(z,3)*exp(y)*sin(M_PI*x) - 1.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x))*cos(MZ*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MX*x)*sin(MY*y) - 0.2*MZ*theta*(-1.5*pow(y,2)*sin(2.0*M_PI*z) + 0.5*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y))*cos(MX*x)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MY*y)*sin(MZ*z) + 2.0*(-2.0*z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 0.5*M_PI*pow(z,3)*cos(M_PI*x)*exp(y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*(z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 2*M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*theta*(1 + 0.1*MY*cos(MX*x)*cos(MY*y)*cos(MZ*z))*(-z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) - 6.0*M_PI*pow(y,2)*cos(2.0*M_PI*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)));
1117:      Fm[2] = -2*z + 8.0*pow(M_PI,2)*(pow(x,3) + pow(y,3))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(2.0*M_PI*z) - 0.2*MX*theta*(-1.5*pow(x,2)*sin(2.0*M_PI*z) + 1.5*pow(z,2)*exp(y)*sin(M_PI*x))*cos(MZ*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MX*x)*sin(MY*y) + 0.4*M_PI*MZ*theta*(pow(x,3) + pow(y,3))*cos(MX*x)*cos(2.0*M_PI*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MY*y)*sin(MZ*z) + 2.0*(-3.0*x*sin(2.0*M_PI*z) + 1.5*M_PI*pow(z,2)*cos(M_PI*x)*exp(y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*(-3.0*y*sin(2.0*M_PI*z) - 0.5*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 0.5*M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*theta*(1 + 0.1*MY*cos(MX*x)*cos(MY*y)*cos(MZ*z))*(-1.5*pow(y,2)*sin(2.0*M_PI*z) + 0.5*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))  ;
1118:      */
1119:     /*
1120:      Fm[0]=-2*x - 2.0*pow(M_PI,2)*sin(M_PI*x);
1121:      Fm[1]=-2*y - 2.0*pow(M_PI,2)*sin(M_PI*y);
1122:      Fm[2]=-2*z - 2.0*pow(M_PI,2)*sin(M_PI*z);
1123:      */
1124:     /*
1125:      Fm[0] = -2*x + pow(z,3)*exp(y)*sin(M_PI*x) + 6.0*z*exp(y)*sin(M_PI*x) - 6.0*M_PI*pow(x,2)*cos(2.0*M_PI*z) - 2.0*pow(M_PI,2)*pow(z,3)*exp(y)*sin(M_PI*x) + 2.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x) - 2.0*z*pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x) ;
1126:      Fm[1] = -2*y - 6.0*M_PI*pow(y,2)*cos(2.0*M_PI*z) + 2.0*z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 6.0*z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*pow(z,3)*cos(M_PI*x)*exp(y) - 4.0*M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y);
1127:      Fm[2] = -2*z - 6.0*x*sin(2.0*M_PI*z) - 6.0*y*sin(2.0*M_PI*z) - cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 8.0*pow(M_PI,2)*(pow(x,3) + pow(y,3))*sin(2.0*M_PI*z) + 3.0*M_PI*pow(z,2)*cos(M_PI*x)*exp(y) + M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;
1128:      */
1129:     Fm[0] = -2*x + 2*z*(PetscPowRealInt(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x) - 1.0*M_PI*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x)) + PetscPowRealInt(z,3)*exp(y)*sin(M_PI*x) + 6.0*z*exp(y)*sin(M_PI*x) - 1.0*PetscPowRealInt(M_PI,2)*PetscPowRealInt(z,3)*exp(y)*sin(M_PI*x) + 2.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x) - 2.0*z*PetscPowRealInt(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x);
1130:     Fm[1] = -2*y + 2*z*(-cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + PetscPowRealInt(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)) + 2.0*z*cos(2.0*M_PI *x)*exp(-y)*sin(M_PI*y) - 6.0*z*PetscPowRealInt(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 4.0*M_PI*z*cos(M_PI *y)*cos(2.0*M_PI*x)*exp(-y);
1131:     Fm[2] = -2*z + PetscPowRealInt(z,2)*(-2.0*PetscPowRealInt(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 2.0*PetscPowRealInt(M_PI,3)*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)) + PetscPowRealInt(z,2)*(cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 - 3*PetscPowRealInt(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + PetscPowRealInt(M_PI,3)*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2 - 3*M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2) + 1.0*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 0.25*PetscPowRealInt(M_PI,3)*PetscPowRealInt(z,4)*cos(M_PI*x)*exp(y) - 0.25*M_PI*PetscPowRealInt(z,4)*cos(M_PI*x)*exp(y) - 3.0*M_PI*PetscPowRealInt(z,2)*cos(M_PI*x)*exp(y) - 1.0*M_PI*cos(M_PI *y)*cos(2.0*M_PI*x)*exp(-y);
1132:   }
1133:   if (Fc) {
1134:     /**Fc = -2.0*M_PI*(pow(x,3) + pow(y,3))*cos(2.0*M_PI*z) - z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*pow(z,3)*cos(M_PI*x)*exp(y) + M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;*/
1135:     /**Fc = M_PI*cos(M_PI*x) + M_PI*cos(M_PI*y) + M_PI*cos(M_PI*z);*/
1136:     /**Fc = -2.0*M_PI*(pow(x,3) + pow(y,3))*cos(2.0*M_PI*z) - z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*pow(z,3)*cos(M_PI*x)*exp(y) + M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y);*/
1137:     *Fc = 0.0;
1138:   }
1139: }

1143: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx,PetscInt my,PetscInt mz,DM *_da,Vec *_X)
1144: {
1145:   DM             da,cda;
1146:   Vec            X,local_X;
1147:   StokesDOF      ***_stokes;
1148:   Vec            coords;
1149:   DMDACoor3d     ***_coords;
1150:   PetscInt       si,sj,sk,ei,ej,ek,i,j,k;

1154:   DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1155:                       mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,4,1,NULL,NULL,NULL,&da);
1156:   DMDASetFieldName(da,0,"anlytic_Vx");
1157:   DMDASetFieldName(da,1,"anlytic_Vy");
1158:   DMDASetFieldName(da,2,"anlytic_Vz");
1159:   DMDASetFieldName(da,3,"analytic_P");

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

1163:   DMGetCoordinatesLocal(da,&coords);
1164:   DMGetCoordinateDM(da,&cda);
1165:   DMDAVecGetArray(cda,coords,&_coords);

1167:   DMCreateGlobalVector(da,&X);
1168:   DMCreateLocalVector(da,&local_X);
1169:   DMDAVecGetArray(da,local_X,&_stokes);

1171:   DMDAGetGhostCorners(da,&si,&sj,&sk,&ei,&ej,&ek);
1172:   for (k = sk; k < sk+ek; k++) {
1173:     for (j = sj; j < sj+ej; j++) {
1174:       for (i = si; i < si+ei; i++) {
1175:         PetscReal pos[NSD],pressure,vel[NSD];

1177:         pos[0] = PetscRealPart(_coords[k][j][i].x);
1178:         pos[1] = PetscRealPart(_coords[k][j][i].y);
1179:         pos[2] = PetscRealPart(_coords[k][j][i].z);

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

1183:         _stokes[k][j][i].u_dof = vel[0];
1184:         _stokes[k][j][i].v_dof = vel[1];
1185:         _stokes[k][j][i].w_dof = vel[2];
1186:         _stokes[k][j][i].p_dof = pressure;
1187:       }
1188:     }
1189:   }
1190:   DMDAVecRestoreArray(da,local_X,&_stokes);
1191:   DMDAVecRestoreArray(cda,coords,&_coords);

1193:   DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1194:   DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);

1196:   VecDestroy(&local_X);

1198:   *_da = da;
1199:   *_X  = X;
1200:   return(0);
1201: }

1205: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da,Vec X,Vec X_analytic)
1206: {
1207:   DM          cda;
1208:   Vec         coords,X_analytic_local,X_local;
1209:   DMDACoor3d  ***_coords;
1210:   PetscInt    sex,sey,sez,mx,my,mz;
1211:   PetscInt    ei,ej,ek;
1212:   PetscScalar el_coords[NODES_PER_EL*NSD];
1213:   StokesDOF   ***stokes_analytic,***stokes;
1214:   StokesDOF   stokes_analytic_e[NODES_PER_EL],stokes_e[NODES_PER_EL];

1216:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1217:   PetscScalar    Ni_p[NODES_PER_EL];
1218:   PetscInt       ngp;
1219:   PetscScalar    gp_xi[GAUSS_POINTS][NSD];
1220:   PetscScalar    gp_weight[GAUSS_POINTS];
1221:   PetscInt       p,i;
1222:   PetscScalar    J_p,fac;
1223:   PetscScalar    h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1224:   PetscScalar    tint_p_ms,tint_p,int_p_ms,int_p;
1225:   PetscInt       M;
1226:   PetscReal      xymin[NSD],xymax[NSD];

1230:   /* define quadrature rule */
1231:   ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

1233:   /* setup for coords */
1234:   DMGetCoordinateDM(stokes_da,&cda);
1235:   DMGetCoordinatesLocal(stokes_da,&coords);
1236:   DMDAVecGetArray(cda,coords,&_coords);

1238:   /* setup for analytic */
1239:   DMCreateLocalVector(stokes_da,&X_analytic_local);
1240:   DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1241:   DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1242:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1244:   /* setup for solution */
1245:   DMCreateLocalVector(stokes_da,&X_local);
1246:   DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1247:   DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1248:   DMDAVecGetArray(stokes_da,X_local,&stokes);

1250:   DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1251:   DMDAGetBoundingBox(stokes_da,xymin,xymax);

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

1255:   tp_L2     = tu_L2 = tu_H1 = 0.0;
1256:   tint_p_ms = tint_p = 0.0;

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

1260:   for (ek = sez; ek < sez+mz; ek++) {
1261:     for (ej = sey; ej < sey+my; ej++) {
1262:       for (ei = sex; ei < sex+mx; ei++) {
1263:         /* get coords for the element */
1264:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1265:         StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1266:         StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);

1268:         /* evaluate integral */
1269:         for (p = 0; p < ngp; p++) {
1270:           ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1271:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1272:           ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1273:           fac = gp_weight[p]*J_p;

1275:           for (i = 0; i < NODES_PER_EL; i++) {
1276:             tint_p_ms = tint_p_ms+fac*Ni_p[i]*stokes_analytic_e[i].p_dof;
1277:             tint_p    = tint_p   +fac*Ni_p[i]*stokes_e[i].p_dof;
1278:           }
1279:         }
1280:       }
1281:     }
1282:   }
1283:   MPI_Allreduce(&tint_p_ms,&int_p_ms,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1284:   MPI_Allreduce(&tint_p,&int_p,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);

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

1288:   /* remove mine and add manufacture one */
1289:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1290:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);

1292:   {
1293:     PetscInt    k,L,dof;
1294:     PetscScalar *fields;
1295:     DMDAGetInfo(stokes_da,0, 0,0,0, 0,0,0, &dof,0,0,0,0,0);

1297:     VecGetLocalSize(X_local,&L);
1298:     VecGetArray(X_local,&fields);
1299:     for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1300:     VecRestoreArray(X_local,&fields);

1302:     VecGetLocalSize(X,&L);
1303:     VecGetArray(X,&fields);
1304:     for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1305:     VecRestoreArray(X,&fields);
1306:   }

1308:   DMDAVecGetArray(stokes_da,X_local,&stokes);
1309:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1311:   for (ek = sez; ek < sez+mz; ek++) {
1312:     for (ej = sey; ej < sey+my; ej++) {
1313:       for (ei = sex; ei < sex+mx; ei++) {
1314:         /* get coords for the element */
1315:         GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1316:         StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1317:         StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);

1319:         /* evaluate integral */
1320:         p_e_L2 = 0.0;
1321:         u_e_L2 = 0.0;
1322:         u_e_H1 = 0.0;
1323:         for (p = 0; p < ngp; p++) {
1324:           ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1325:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1326:           ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1327:           fac = gp_weight[p]*J_p;

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

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

1334:             u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1335:             v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1336:             w_error = stokes_e[i].w_dof-stokes_analytic_e[i].w_dof;
1337:             /*
1338:              if (p==0) {
1339:              printf("p=0: %d %d %d %1.4e,%1.4e,%1.4e \n", ei,ej,ek,u_error,v_error,w_error);
1340:              }
1341:              */
1342:             u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error+w_error*w_error);

1344:             u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error              /* du/dx */
1345:                                  +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error
1346:                                  +GNx_p[2][i]*u_error*GNx_p[2][i]*u_error
1347:                                  +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error               /* dv/dx */
1348:                                  +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error
1349:                                  +GNx_p[2][i]*v_error*GNx_p[2][i]*v_error
1350:                                  +GNx_p[0][i]*w_error*GNx_p[0][i]*w_error               /* dw/dx */
1351:                                  +GNx_p[1][i]*w_error*GNx_p[1][i]*w_error
1352:                                  +GNx_p[2][i]*w_error*GNx_p[2][i]*w_error);
1353:           }
1354:         }

1356:         tp_L2 += p_e_L2;
1357:         tu_L2 += u_e_L2;
1358:         tu_H1 += u_e_H1;
1359:       }
1360:     }
1361:   }
1362:   MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1363:   MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1364:   MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1365:   p_L2 = PetscSqrtScalar(p_L2);
1366:   u_L2 = PetscSqrtScalar(u_L2);
1367:   u_H1 = PetscSqrtScalar(u_H1);

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


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

1374:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1375:   VecDestroy(&X_analytic_local);
1376:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1377:   VecDestroy(&X_local);
1378:   return(0);
1379: }

1383: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da,Vec FIELD,const char file_prefix[])
1384: {
1385:   char           vtk_filename[PETSC_MAX_PATH_LEN];
1386:   PetscMPIInt    rank;
1387:   MPI_Comm       comm;
1388:   FILE           *vtk_fp = NULL;
1389:   PetscInt       si,sj,sk,nx,ny,nz,i;
1390:   PetscInt       f,n_fields,N;
1391:   DM             cda;
1392:   Vec            coords;
1393:   Vec            l_FIELD;
1394:   PetscScalar    *_L_FIELD;
1395:   PetscInt       memory_offset;
1396:   PetscScalar    *buffer;
1397:   PetscLogDouble t0,t1;

1401:   PetscTime(&t0);

1403:   /* create file name */
1404:   PetscObjectGetComm((PetscObject)da,&comm);
1405:   MPI_Comm_rank(comm,&rank);
1406:   PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"subdomain-%s-p%1.4d.vts",file_prefix,rank);

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

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

1414:   /* coords */
1415:   DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1416:   N    = nx * ny * nz;

1418: #if defined(PETSC_WORDS_BIGENDIAN)
1419:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1420: #else
1421:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1422: #endif
1423:   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);
1424:   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);

1426:   memory_offset = 0;

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

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

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

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

1440:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      <PointData Scalars=\" ");
1441:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_fields,0,0,0,0,0);
1442:   for (f=0; f<n_fields; f++) {
1443:     const char *field_name;
1444:     DMDAGetFieldName(da,f,&field_name);
1445:     PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"%s ",field_name);
1446:   }
1447:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\">\n");

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

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

1457:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      </PointData>\n");
1458:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </Piece>\n");
1459:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  </StructuredGrid>\n");

1461:   PetscMalloc(sizeof(PetscScalar)*N,&buffer);
1462:   DMGetLocalVector(da,&l_FIELD);
1463:   DMGlobalToLocalBegin(da, FIELD,INSERT_VALUES,l_FIELD);
1464:   DMGlobalToLocalEnd(da,FIELD,INSERT_VALUES,l_FIELD);
1465:   VecGetArray(l_FIELD,&_L_FIELD);

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

1470:   /* write coordinates */
1471:   {
1472:     int         length = sizeof(PetscScalar)*N*3;
1473:     PetscScalar *allcoords;

1475:     fwrite(&length,sizeof(int),1,vtk_fp);
1476:     VecGetArray(coords,&allcoords);
1477:     fwrite(allcoords,sizeof(PetscScalar),3*N,vtk_fp);
1478:     VecRestoreArray(coords,&allcoords);
1479:   }
1480:   /* write fields */
1481:   for (f=0; f<n_fields; f++) {
1482:     int length = sizeof(PetscScalar)*N;
1483:     fwrite(&length,sizeof(int),1,vtk_fp);
1484:     /* load */
1485:     for (i=0; i<N; i++) buffer[i] = _L_FIELD[n_fields*i + f];

1487:     /* write */
1488:     fwrite(buffer,sizeof(PetscScalar),N,vtk_fp);
1489:   }
1490:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\n  </AppendedData>\n");

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

1494:   PetscFree(buffer);
1495:   VecRestoreArray(l_FIELD,&_L_FIELD);
1496:   DMRestoreLocalVector(da,&l_FIELD);

1498:   if (vtk_fp) {
1499:     fclose(vtk_fp);
1500:     vtk_fp = NULL;
1501:   }

1503:   PetscTime(&t1);
1504:   return(0);
1505: }

1509: PetscErrorCode DAViewVTK_write_PieceExtend(FILE *vtk_fp,PetscInt indent_level,DM da,const char local_file_prefix[])
1510: {
1511:   PetscMPIInt    nproc,rank;
1512:   MPI_Comm       comm;
1513:   const PetscInt *lx,*ly,*lz;
1514:   PetscInt       M,N,P,pM,pN,pP,sum,*olx,*oly,*olz;
1515:   PetscInt       *osx,*osy,*osz,*oex,*oey,*oez;
1516:   PetscInt       i,j,k,II,stencil;

1520:   /* create file name */
1521:   PetscObjectGetComm((PetscObject)da,&comm);
1522:   MPI_Comm_size(comm,&nproc);
1523:   MPI_Comm_rank(comm,&rank);

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

1528:   /* generate start,end list */
1529:   PetscMalloc(sizeof(PetscInt)*(pM+1),&olx);
1530:   PetscMalloc(sizeof(PetscInt)*(pN+1),&oly);
1531:   PetscMalloc(sizeof(PetscInt)*(pP+1),&olz);
1532:   sum  = 0;
1533:   for (i=0; i<pM; i++) {
1534:     olx[i] = sum;
1535:     sum    = sum + lx[i];
1536:   }
1537:   olx[pM] = sum;
1538:   sum     = 0;
1539:   for (i=0; i<pN; i++) {
1540:     oly[i] = sum;
1541:     sum    = sum + ly[i];
1542:   }
1543:   oly[pN] = sum;
1544:   sum     = 0;
1545:   for (i=0; i<pP; i++) {
1546:     olz[i] = sum;
1547:     sum    = sum + lz[i];
1548:   }
1549:   olz[pP] = sum;

1551:   PetscMalloc(sizeof(PetscInt)*(pM),&osx);
1552:   PetscMalloc(sizeof(PetscInt)*(pN),&osy);
1553:   PetscMalloc(sizeof(PetscInt)*(pP),&osz);
1554:   PetscMalloc(sizeof(PetscInt)*(pM),&oex);
1555:   PetscMalloc(sizeof(PetscInt)*(pN),&oey);
1556:   PetscMalloc(sizeof(PetscInt)*(pP),&oez);
1557:   for (i=0; i<pM; i++) {
1558:     osx[i] = olx[i] - stencil;
1559:     oex[i] = olx[i] + lx[i] + stencil;
1560:     if (osx[i]<0) osx[i]=0;
1561:     if (oex[i]>M) oex[i]=M;
1562:   }

1564:   for (i=0; i<pN; i++) {
1565:     osy[i] = oly[i] - stencil;
1566:     oey[i] = oly[i] + ly[i] + stencil;
1567:     if (osy[i]<0)osy[i]=0;
1568:     if (oey[i]>M)oey[i]=N;
1569:   }
1570:   for (i=0; i<pP; i++) {
1571:     osz[i] = olz[i] - stencil;
1572:     oez[i] = olz[i] + lz[i] + stencil;
1573:     if (osz[i]<0) osz[i]=0;
1574:     if (oez[i]>P) oez[i]=P;
1575:   }

1577:   for (k=0; k<pP; k++) {
1578:     for (j=0; j<pN; j++) {
1579:       for (i=0; i<pM; i++) {
1580:         char     name[PETSC_MAX_PATH_LEN];
1581:         PetscInt procid = i + j*pM + k*pM*pN; /* convert proc(i,j,k) to pid */
1582:         PetscSNPrintf(name,sizeof(name),"subdomain-%s-p%1.4d.vts",local_file_prefix,procid);
1583:         for (II=0; II<indent_level; II++) PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  ");

1585:         PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<Piece Extent=\"%d %d %d %d %d %d\"      Source=\"%s\"/>\n",
1586:                      osx[i],oex[i]-1,
1587:                      osy[j],oey[j]-1,
1588:                      osz[k],oez[k]-1,name);
1589:       }
1590:     }
1591:   }
1592:   PetscFree(olx);
1593:   PetscFree(oly);
1594:   PetscFree(olz);
1595:   PetscFree(osx);
1596:   PetscFree(osy);
1597:   PetscFree(osz);
1598:   PetscFree(oex);
1599:   PetscFree(oey);
1600:   PetscFree(oez);
1601:   return(0);
1602: }

1606: PetscErrorCode DAView_3DVTK_PStructuredGrid(DM da,const char file_prefix[],const char local_file_prefix[])
1607: {
1608:   MPI_Comm       comm;
1609:   PetscMPIInt    nproc,rank;
1610:   char           vtk_filename[PETSC_MAX_PATH_LEN];
1611:   FILE           *vtk_fp = NULL;
1612:   PetscInt       M,N,P,si,sj,sk,nx,ny,nz;
1613:   PetscInt       i,dofs;

1617:   /* only master generates this file */
1618:   PetscObjectGetComm((PetscObject)da,&comm);
1619:   MPI_Comm_size(comm,&nproc);
1620:   MPI_Comm_rank(comm,&rank);

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

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

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

1632: #if defined(PETSC_WORDS_BIGENDIAN)
1633:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1634: #else
1635:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1636: #endif

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

1643:   /* DUMP THE CELL REFERENCES */
1644:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    <PCellData>\n");
1645:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </PCellData>\n");

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

1651:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    <PPointData>\n");
1652:   for (i=0; i<dofs; i++) {
1653:     const char *fieldname;
1654:     DMDAGetFieldName(da,i,&fieldname);
1655:     PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"      <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n",fieldname);
1656:   }
1657:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"    </PPointData>\n");

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

1662:   /* close the file */
1663:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"  </PStructuredGrid>\n");
1664:   PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");

1666:   if (vtk_fp) {
1667:     fclose(vtk_fp);
1668:     vtk_fp = NULL;
1669:   }
1670:   return(0);
1671: }

1675: PetscErrorCode DAView3DPVTS(DM da, Vec x,const char NAME[])
1676: {
1677:   char           vts_filename[PETSC_MAX_PATH_LEN];
1678:   char           pvts_filename[PETSC_MAX_PATH_LEN];

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

1685:   PetscSNPrintf(pvts_filename,sizeof(pvts_filename),"%s-mesh",NAME);
1686:   DAView_3DVTK_PStructuredGrid(da,pvts_filename,vts_filename);
1687:   return(0);
1688: }

1692: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
1693: {
1695:   PetscReal      norms[4];
1696:   Vec            Br,v,w;
1697:   Mat            A;

1700:   KSPGetOperators(ksp,&A,0,0);
1701:   MatGetVecs(A,&w,&v);

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

1705:   VecStrideNorm(Br,0,NORM_2,&norms[0]);
1706:   VecStrideNorm(Br,1,NORM_2,&norms[1]);
1707:   VecStrideNorm(Br,2,NORM_2,&norms[2]);
1708:   VecStrideNorm(Br,3,NORM_2,&norms[3]);

1710:   VecDestroy(&v);
1711:   VecDestroy(&w);

1713:   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]);
1714:   return(0);
1715: }

1719: static PetscErrorCode PCMGSetupViaCoarsen(PC pc,DM da_fine)
1720: {
1721:   PetscInt       nlevels,k,PETSC_UNUSED finest;
1722:   DM             *da_list,*daclist;
1723:   Mat            R;

1727:   nlevels = 1;
1728:   PetscOptionsGetInt(NULL,"-levels",&nlevels,0);

1730:   PetscMalloc(sizeof(DM)*nlevels,&da_list);
1731:   for (k=0; k<nlevels; k++) da_list[k] = NULL;
1732:   PetscMalloc(sizeof(DM)*nlevels,&daclist);
1733:   for (k=0; k<nlevels; k++) daclist[k] = NULL;

1735:   /* finest grid is nlevels - 1 */
1736:   finest     = nlevels - 1;
1737:   daclist[0] = da_fine;
1738:   PetscObjectReference((PetscObject)da_fine);
1739:   DMCoarsenHierarchy(da_fine,nlevels-1,&daclist[1]);
1740:   for (k=0; k<nlevels; k++) {
1741:     da_list[k] = daclist[nlevels-1-k];
1742:     DMDASetUniformCoordinates(da_list[k],0.0,1.0,0.0,1.0,0.0,1.0);
1743:   }

1745:   PCMGSetLevels(pc,nlevels,NULL);
1746:   PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
1747:   PCMGSetGalerkin(pc,PETSC_TRUE);

1749:   for (k=1; k<nlevels; k++) {
1750:     DMCreateInterpolation(da_list[k-1],da_list[k],&R,NULL);
1751:     PCMGSetInterpolation(pc,k,R);
1752:     MatDestroy(&R);
1753:   }

1755:   /* tidy up */
1756:   for (k=0; k<nlevels; k++) {
1757:     DMDestroy(&da_list[k]);
1758:   }
1759:   PetscFree(da_list);
1760:   PetscFree(daclist);
1761:   return(0);
1762: }

1766: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx,PetscInt my,PetscInt mz)
1767: {
1768:   DM             da_Stokes;
1769:   PetscInt       u_dof,p_dof,dof,stencil_width;
1770:   Mat            A,B;
1771:   PetscInt       mxl,myl,mzl;
1772:   DM             vel_cda;
1773:   Vec            vel_coords;
1774:   PetscInt       p;
1775:   Vec            f,X;
1776:   DMDACoor3d     ***_vel_coords;
1777:   PetscInt       its;
1778:   KSP            ksp_S;
1779:   PetscInt       model_definition = 0;
1780:   PetscInt       ei,ej,ek,sex,sey,sez,Imx,Jmy,Kmz;
1781:   CellProperties cell_properties;
1782:   PetscBool      write_output = PETSC_FALSE;

1786:   /* Generate the da for velocity and pressure */
1787:   /* Num nodes in each direction is mx+1, my+1, mz+1 */
1788:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1789:   p_dof         = P_DOFS; /* p - pressure */
1790:   dof           = u_dof+p_dof;
1791:   stencil_width = 1;
1792:   DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1793:                                mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,NULL,&da_Stokes);
1794:   DMDASetFieldName(da_Stokes,0,"Vx");
1795:   DMDASetFieldName(da_Stokes,1,"Vy");
1796:   DMDASetFieldName(da_Stokes,2,"Vz");
1797:   DMDASetFieldName(da_Stokes,3,"P");

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

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

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

1808:   /* interpolate the coordinates to quadrature points */
1809:   DMGetCoordinateDM(da_Stokes,&vel_cda);
1810:   DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1811:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1812:   DMDAGetElementCorners(da_Stokes,&sex,&sey,&sez,&Imx,&Jmy,&Kmz);
1813:   for (ek = sez; ek < sez+Kmz; ek++) {
1814:     for (ej = sey; ej < sey+Jmy; ej++) {
1815:       for (ei = sex; ei < sex+Imx; ei++) {
1816:         /* get coords for the element */
1817:         PetscInt               ngp;
1818:         PetscScalar            gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
1819:         PetscScalar            el_coords[NSD*NODES_PER_EL];
1820:         GaussPointCoefficients *cell;

1822:         CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1823:         GetElementCoords3D(_vel_coords,ei,ej,ek,el_coords);
1824:         ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);

1826:         for (p = 0; p < GAUSS_POINTS; p++) {
1827:           PetscScalar xi_p[NSD],Ni_p[NODES_PER_EL];
1828:           PetscScalar gp_x,gp_y,gp_z;
1829:           PetscInt    n;

1831:           xi_p[0] = gp_xi[p][0];
1832:           xi_p[1] = gp_xi[p][1];
1833:           xi_p[2] = gp_xi[p][2];
1834:           ShapeFunctionQ13D_Evaluate(xi_p,Ni_p);

1836:           gp_x = gp_y = gp_z = 0.0;
1837:           for (n = 0; n < NODES_PER_EL; n++) {
1838:             gp_x = gp_x+Ni_p[n]*el_coords[NSD*n];
1839:             gp_y = gp_y+Ni_p[n]*el_coords[NSD*n+1];
1840:             gp_z = gp_z+Ni_p[n]*el_coords[NSD*n+2];
1841:           }
1842:           cell->gp_coords[NSD*p]   = gp_x;
1843:           cell->gp_coords[NSD*p+1] = gp_y;
1844:           cell->gp_coords[NSD*p+2] = gp_z;
1845:         }
1846:       }
1847:     }
1848:   }

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

1852:   switch (model_definition) {
1853:   case 0: /* isoviscous */
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:           GaussPointCoefficients *cell;

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

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

1867:             cell->fx[p] = 0.0*coord_x;
1868:             cell->fy[p] = -sin((double)2.2*M_PI*coord_y)*cos(1.0*M_PI*coord_x);
1869:             cell->fz[p] = 0.0*coord_z;
1870:             cell->hc[p] = 0.0;
1871:           }
1872:         }
1873:       }
1874:     }
1875:     break;

1877:   case 1: /* manufactured */
1878:     for (ek = sez; ek < sez+Kmz; ek++) {
1879:       for (ej = sey; ej < sey+Jmy; ej++) {
1880:         for (ei = sex; ei < sex+Imx; ei++) {
1881:           PetscReal              eta,Fm[NSD],Fc,pos[NSD];
1882:           GaussPointCoefficients *cell;

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

1890:             pos[0] = coord_x;
1891:             pos[1] = coord_y;
1892:             pos[2] = coord_z;

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

1897:             cell->fx[p] = Fm[0];
1898:             cell->fy[p] = Fm[1];
1899:             cell->fz[p] = Fm[2];
1900:             cell->hc[p] = Fc;
1901:           }
1902:         }
1903:       }
1904:     }
1905:     break;

1907:   case 2: /* solcx */
1908:     for (ek = sez; ek < sez+Kmz; ek++) {
1909:       for (ej = sey; ej < sey+Jmy; ej++) {
1910:         for (ei = sex; ei < sex+Imx; ei++) {
1911:           GaussPointCoefficients *cell;

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

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

1921:             cell->fx[p] = 0.0;
1922:             cell->fy[p] = -sin((double)3*M_PI*coord_y)*cos(1.0*M_PI*coord_x);
1923:             cell->fz[p] = 0.0*coord_z;
1924:             cell->hc[p] = 0.0;
1925:           }
1926:         }
1927:       }
1928:     }
1929:     break;

1931:   case 3: /* sinker */
1932:     for (ek = sez; ek < sez+Kmz; ek++) {
1933:       for (ej = sey; ej < sey+Jmy; ej++) {
1934:         for (ei = sex; ei < sex+Imx; ei++) {
1935:           GaussPointCoefficients *cell;

1937:           CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1938:           for (p = 0; p < GAUSS_POINTS; p++) {
1939:             PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p]);
1940:             PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);
1941:             PetscReal zp = PetscRealPart(cell->gp_coords[NSD*p+2]);

1943:             cell->eta[p] = 1.0e-2;
1944:             cell->fx[p]  = 0.0;
1945:             cell->fy[p]  = 0.0;
1946:             cell->fz[p]  = 0.0;
1947:             cell->hc[p]  = 0.0;

1949:             if ((fabs(xp-0.5) < 0.2) &&
1950:                 (fabs(yp-0.5) < 0.2) &&
1951:                 (fabs(zp-0.5) < 0.2)) {
1952:               cell->eta[p] = 1.0;
1953:               cell->fz[p]  = 1.0;
1954:             }

1956:           }
1957:         }
1958:       }
1959:     }
1960:     break;

1962:   default:
1963:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"No default model is supported. Choose either -model {0,1,2,3}");
1964:     break;
1965:   }

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

1969:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1970:   DMCreateMatrix(da_Stokes,MATAIJ,&A);
1971:   DMCreateMatrix(da_Stokes,MATAIJ,&B);
1972:   DMCreateGlobalVector(da_Stokes,&X);
1973:   DMCreateGlobalVector(da_Stokes,&f);

1975:   /* assemble A11 */
1976:   MatZeroEntries(A);
1977:   MatZeroEntries(B);
1978:   VecZeroEntries(f);

1980:   AssembleA_Stokes(A,da_Stokes,cell_properties);
1981:   AssembleA_PCStokes(B,da_Stokes,cell_properties);
1982:   /* build force vector */
1983:   AssembleF_Stokes(f,da_Stokes,cell_properties);

1985:   /* SOLVE */
1986:   KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1987:   KSPSetOptionsPrefix(ksp_S,"stokes_"); /* stokes */
1988:   KSPSetOperators(ksp_S,A,B,SAME_NONZERO_PATTERN);
1989:   KSPSetFromOptions(ksp_S);

1991:   {
1992:     PC             pc;
1993:     const PetscInt ufields[] = {0,1,2},pfields[] = {3};
1994:     KSPGetPC(ksp_S,&pc);
1995:     PCFieldSplitSetBlockSize(pc,4);
1996:     PCFieldSplitSetFields(pc,"u",3,ufields,ufields);
1997:     PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1998:   }

2000:   {
2001:     PC        pc;
2002:     PetscBool same = PETSC_FALSE;
2003:     KSPGetPC(ksp_S,&pc);
2004:     PetscObjectTypeCompare((PetscObject)pc,PCMG,&same);
2005:     if (same) {
2006:       PCMGSetupViaCoarsen(pc,da_Stokes);
2007:     }
2008:   }

2010:   {
2011:     PetscBool stokes_monitor = PETSC_FALSE;
2012:     PetscOptionsGetBool(NULL,"-stokes_ksp_monitor_blocks",&stokes_monitor,0);
2013:     if (stokes_monitor) {
2014:       KSPMonitorSet(ksp_S,KSPMonitorStokesBlocks,NULL,NULL);
2015:     }
2016:   }
2017:   KSPSolve(ksp_S,f,X);

2019:   PetscOptionsGetBool(NULL,"-write_pvts",&write_output,NULL);
2020:   if (write_output) {DAView3DPVTS(da_Stokes,X,"up");}
2021:   {
2022:     PetscBool flg = PETSC_FALSE;
2023:     char      filename[PETSC_MAX_PATH_LEN];
2024:     PetscOptionsGetString(NULL,"-write_binary",filename,sizeof(filename),&flg);
2025:     if (flg) {
2026:       PetscViewer viewer;
2027:       /* PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer); */
2028:       PetscViewerVTKOpen(PETSC_COMM_WORLD,"ex42.vts",FILE_MODE_WRITE,&viewer);
2029:       VecView(X,viewer);
2030:       PetscViewerDestroy(&viewer);
2031:     }
2032:   }
2033:   KSPGetIterationNumber(ksp_S,&its);

2035:   /* verify */
2036:   if (model_definition == 1) {
2037:     DM  da_Stokes_analytic;
2038:     Vec X_analytic;

2040:     DMDACreateManufacturedSolution(mx,my,mz,&da_Stokes_analytic,&X_analytic);
2041:     if (write_output) {
2042:       DAView3DPVTS(da_Stokes_analytic,X_analytic,"ms");
2043:     }
2044:     DMDAIntegrateErrors3D(da_Stokes_analytic,X,X_analytic);
2045:     if (write_output) {
2046:       DAView3DPVTS(da_Stokes,X,"up2");
2047:     }
2048:     DMDestroy(&da_Stokes_analytic);
2049:     VecDestroy(&X_analytic);
2050:   }

2052:   KSPDestroy(&ksp_S);
2053:   VecDestroy(&X);
2054:   VecDestroy(&f);
2055:   MatDestroy(&A);
2056:   MatDestroy(&B);

2058:   CellPropertiesDestroy(&cell_properties);
2059:   DMDestroy(&da_Stokes);
2060:   return(0);
2061: }

2065: int main(int argc,char **args)
2066: {
2068:   PetscInt       mx,my,mz;

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

2072:   mx   = my = mz = 10;
2073:   PetscOptionsGetInt(NULL,"-mx",&mx,NULL);
2074:   my   = mx; mz = mx;
2075:   PetscOptionsGetInt(NULL,"-my",&my,NULL);
2076:   PetscOptionsGetInt(NULL,"-mz",&mz,NULL);

2078:   solve_stokes_3d_coupled(mx,my,mz);

2080:   PetscFinalize();
2081:   return 0;
2082: }