Actual source code: ex42.c
1: static char help[] =
2: "Solves the incompressible, variable viscosity stokes equation in 3d using Q1Q1 elements, \n\
3: stabilized with Bochev's polynomial projection method. Note that implementation here assumes \n\
4: all boundaries are free-slip, i.e. zero normal flow and zero tangential stress \n\
5: -mx : number elements in x-direction \n\
6: -my : number elements in y-direction \n\
7: -mz : number elements in z-direction \n\
8: -stokes_ksp_monitor_blocks : active monitor for each component u,v,w,p \n\
9: -model : defines viscosity and forcing function. Choose either: 0 (isoviscous), 1 (manufactured-broken), 2 (solcx-2d), 3 (sinker) \n";
11: /* Contributed by Dave May */
13: #include petscksp.h
14: #include petscpcmg.h
15: #include petscdmda.h
17: #define PROFILE_TIMING
18: #define ASSEMBLE_LOWER_TRIANGULAR
20: #define NSD 3 /* number of spatial dimensions */
21: #define NODES_PER_EL 8 /* nodes per element */
22: #define U_DOFS 3 /* degrees of freedom per velocity node */
23: #define P_DOFS 1 /* degrees of freedom per pressure node */
24: #define GAUSS_POINTS 8
26: /* Gauss point based evaluation */
27: typedef struct {
28: PetscScalar gp_coords[NSD*GAUSS_POINTS];
29: PetscScalar eta[GAUSS_POINTS];
30: PetscScalar fx[GAUSS_POINTS];
31: PetscScalar fy[GAUSS_POINTS];
32: PetscScalar fz[GAUSS_POINTS];
33: PetscScalar hc[GAUSS_POINTS];
34: } GaussPointCoefficients;
36: typedef struct {
37: PetscScalar u_dof;
38: PetscScalar v_dof;
39: PetscScalar w_dof;
40: PetscScalar p_dof;
41: } StokesDOF;
44: /* FEM routines */
45: /*
46: Element: Local basis function ordering
47: 1-----2
48: | |
49: | |
50: 0-----3
51: */
52: static void ShapeFunctionQ13D_Evaluate(PetscScalar _xi[],PetscScalar Ni[])
53: {
54: PetscReal xi = PetscRealPart(_xi[0]);
55: PetscReal eta = PetscRealPart(_xi[1]);
56: PetscReal zeta = PetscRealPart(_xi[2]);
58: Ni[0] = 0.125 * ( 1.0 - xi ) * ( 1.0 - eta ) * ( 1.0 - zeta );
59: Ni[1] = 0.125 * ( 1.0 - xi ) * ( 1.0 + eta ) * ( 1.0 - zeta );
60: Ni[2] = 0.125 * ( 1.0 + xi ) * ( 1.0 + eta ) * ( 1.0 - zeta );
61: Ni[3] = 0.125 * ( 1.0 + xi ) * ( 1.0 - eta ) * ( 1.0 - zeta );
63: Ni[4] = 0.125 * ( 1.0 - xi ) * ( 1.0 - eta ) * ( 1.0 + zeta );
64: Ni[5] = 0.125 * ( 1.0 - xi ) * ( 1.0 + eta ) * ( 1.0 + zeta );
65: Ni[6] = 0.125 * ( 1.0 + xi ) * ( 1.0 + eta ) * ( 1.0 + zeta );
66: Ni[7] = 0.125 * ( 1.0 + xi ) * ( 1.0 - eta ) * ( 1.0 + zeta );
67: }
69: static void ShapeFunctionQ13D_Evaluate_dxi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
70: {
71: PetscReal xi = PetscRealPart(_xi[0]);
72: PetscReal eta = PetscRealPart(_xi[1]);
73: PetscReal zeta = PetscRealPart(_xi[2]);
74: /* xi deriv */
75: GNi[0][0] = - 0.125 * ( 1.0 - eta ) * ( 1.0 - zeta );
76: GNi[0][1] = - 0.125 * ( 1.0 + eta ) * ( 1.0 - zeta );
77: GNi[0][2] = 0.125 * ( 1.0 + eta ) * ( 1.0 - zeta );
78: GNi[0][3] = 0.125 * ( 1.0 - eta ) * ( 1.0 - zeta );
80: GNi[0][4] = - 0.125 * ( 1.0 - eta ) * ( 1.0 + zeta );
81: GNi[0][5] = - 0.125 * ( 1.0 + eta ) * ( 1.0 + zeta );
82: GNi[0][6] = 0.125 * ( 1.0 + eta ) * ( 1.0 + zeta );
83: GNi[0][7] = 0.125 * ( 1.0 - eta ) * ( 1.0 + zeta );
84: /* eta deriv */
85: GNi[1][0] = - 0.125 * ( 1.0 - xi ) * ( 1.0 - zeta );
86: GNi[1][1] = 0.125 * ( 1.0 - xi ) * ( 1.0 - zeta );
87: GNi[1][2] = 0.125 * ( 1.0 + xi ) * ( 1.0 - zeta );
88: GNi[1][3] = - 0.125 * ( 1.0 + xi ) * ( 1.0 - zeta );
90: GNi[1][4] = - 0.125 * ( 1.0 - xi ) * ( 1.0 + zeta );
91: GNi[1][5] = 0.125 * ( 1.0 - xi ) * ( 1.0 + zeta );
92: GNi[1][6] = 0.125 * ( 1.0 + xi ) * ( 1.0 + zeta );
93: GNi[1][7] = - 0.125 * ( 1.0 + xi ) * ( 1.0 + zeta );
94: /* zeta deriv */
95: GNi[2][0] = -0.125 * ( 1.0 - xi ) * ( 1.0 - eta );
96: GNi[2][1] = -0.125 * ( 1.0 - xi ) * ( 1.0 + eta );
97: GNi[2][2] = -0.125 * ( 1.0 + xi ) * ( 1.0 + eta );
98: GNi[2][3] = -0.125 * ( 1.0 + xi ) * ( 1.0 - eta );
100: GNi[2][4] = 0.125 * ( 1.0 - xi ) * ( 1.0 - eta );
101: GNi[2][5] = 0.125 * ( 1.0 - xi ) * ( 1.0 + eta );
102: GNi[2][6] = 0.125 * ( 1.0 + xi ) * ( 1.0 + eta );
103: GNi[2][7] = 0.125 * ( 1.0 + xi ) * ( 1.0 - eta );
104: }
106: static void matrix_inverse_3x3(PetscScalar A[3][3],PetscScalar B[3][3])
107: {
108: PetscScalar t4, t6, t8, t10, t12, t14, t17;
110: t4 = A[2][0] * A[0][1];
111: t6 = A[2][0] * A[0][2];
112: t8 = A[1][0] * A[0][1];
113: t10 = A[1][0] * A[0][2];
114: t12 = A[0][0] * A[1][1];
115: t14 = A[0][0] * A[1][2];
116: 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]);
118: B[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) * t17;
119: B[0][1] = -(A[0][1] * A[2][2] - A[0][2] * A[2][1]) * t17;
120: B[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) * t17;
121: B[1][0] = -(-A[2][0] * A[1][2] + A[1][0] * A[2][2]) * t17;
122: B[1][1] = (-t6 + A[0][0] * A[2][2]) * t17;
123: B[1][2] = -(-t10 + t14) * t17;
124: B[2][0] = (-A[2][0] * A[1][1] + A[1][0] * A[2][1]) * t17;
125: B[2][1] = -(-t4 + A[0][0] * A[2][1]) * t17;
126: B[2][2] = (-t8 + t12) * t17;
127: }
129: static void ShapeFunctionQ13D_Evaluate_dx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
130: {
131: PetscScalar J00,J01,J02,J10,J11,J12,J20,J21,J22;
132: PetscInt n;
133: PetscScalar iJ[3][3],JJ[3][3];
135: J00 = J01 = J02 = 0.0;
136: J10 = J11 = J12 = 0.0;
137: J20 = J21 = J22 = 0.0;
138: for (n=0; n<NODES_PER_EL; n++) {
139: PetscScalar cx = coords[ NSD*n + 0 ];
140: PetscScalar cy = coords[ NSD*n + 1 ];
141: PetscScalar cz = coords[ NSD*n + 2 ];
143: /* J_ij = d(x_j) / d(xi_i) *//* J_ij = \sum _I GNi[j][I} * x_i */
144: J00 = J00 + GNi[0][n] * cx; /* J_xx */
145: J01 = J01 + GNi[0][n] * cy; /* J_xy = dx/deta */
146: J02 = J02 + GNi[0][n] * cz; /* J_xz = dx/dzeta */
148: J10 = J10 + GNi[1][n] * cx; /* J_yx = dy/dxi */
149: J11 = J11 + GNi[1][n] * cy; /* J_yy */
150: J12 = J12 + GNi[1][n] * cz; /* J_yz */
152: J20 = J20 + GNi[2][n] * cx; /* J_zx */
153: J21 = J21 + GNi[2][n] * cy; /* J_zy */
154: J22 = J22 + GNi[2][n] * cz; /* J_zz */
155: }
157: JJ[0][0] = J00; JJ[0][1] = J01; JJ[0][2] = J02;
158: JJ[1][0] = J10; JJ[1][1] = J11; JJ[1][2] = J12;
159: JJ[2][0] = J20; JJ[2][1] = J21; JJ[2][2] = J22;
161: matrix_inverse_3x3(JJ,iJ);
163: *det_J = J00*J11*J22 - J00*J12*J21 - J10*J01*J22 + J10*J02*J21 + J20*J01*J12 - J20*J02*J11;
165: for (n=0; n<NODES_PER_EL; n++) {
166: GNx[0][n] = GNi[0][n]*iJ[0][0] + GNi[1][n]*iJ[0][1] + GNi[2][n]*iJ[0][2];
167: GNx[1][n] = GNi[0][n]*iJ[1][0] + GNi[1][n]*iJ[1][1] + GNi[2][n]*iJ[1][2];
168: GNx[2][n] = GNi[0][n]*iJ[2][0] + GNi[1][n]*iJ[2][1] + GNi[2][n]*iJ[2][2];
169: }
170: }
172: static void ConstructGaussQuadrature3D(PetscInt *ngp,PetscScalar gp_xi[][NSD],PetscScalar gp_weight[])
173: {
174: *ngp = 8;
175: gp_xi[0][0] = -0.57735026919; gp_xi[0][1] = -0.57735026919; gp_xi[0][2] = -0.57735026919;
176: gp_xi[1][0] = -0.57735026919; gp_xi[1][1] = 0.57735026919; gp_xi[1][2] = -0.57735026919;
177: gp_xi[2][0] = 0.57735026919; gp_xi[2][1] = 0.57735026919; gp_xi[2][2] = -0.57735026919;
178: gp_xi[3][0] = 0.57735026919; gp_xi[3][1] = -0.57735026919; gp_xi[3][2] = -0.57735026919;
180: gp_xi[4][0] = -0.57735026919; gp_xi[4][1] = -0.57735026919; gp_xi[4][2] = 0.57735026919;
181: gp_xi[5][0] = -0.57735026919; gp_xi[5][1] = 0.57735026919; gp_xi[5][2] = 0.57735026919;
182: gp_xi[6][0] = 0.57735026919; gp_xi[6][1] = 0.57735026919; gp_xi[6][2] = 0.57735026919;
183: gp_xi[7][0] = 0.57735026919; gp_xi[7][1] = -0.57735026919; gp_xi[7][2] = 0.57735026919;
185: gp_weight[0] = 1.0;
186: gp_weight[1] = 1.0;
187: gp_weight[2] = 1.0;
188: gp_weight[3] = 1.0;
190: gp_weight[4] = 1.0;
191: gp_weight[5] = 1.0;
192: gp_weight[6] = 1.0;
193: gp_weight[7] = 1.0;
194: }
196: /* procs to the left claim the ghost node as their element */
199: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
200: {
201: PetscInt m,n,p,M,N,P;
202: PetscInt sx,sy,sz;
206: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
207: DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);
209: if (mxl) {
210: *mxl = m;
211: if ((sx+m) == M) { /* last proc */
212: *mxl = m-1;
213: }
214: }
215: if (myl) {
216: *myl = n;
217: if ((sy+n) == N) { /* last proc */
218: *myl = n-1;
219: }
220: }
221: if (mzl) {
222: *mzl = p;
223: if ((sz+p) == P) { /* last proc */
224: *mzl = p-1;
225: }
226: }
227: return(0);
228: }
232: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
233: {
234: PetscInt si,sj,sk;
238: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
240: if (sx) {
241: *sx = si;
242: if (si != 0) {
243: *sx = si+1;
244: }
245: }
246: if (sy) {
247: *sy = sj;
248: if (sj != 0) {
249: *sy = sj+1;
250: }
251: }
252: if (sz) {
253: *sz = sk;
254: if (sk != 0) {
255: *sz = sk+1;
256: }
257: }
258: DMDAGetLocalElementSize(da,mx,my,mz);
259: return(0);
260: }
264: static PetscErrorCode DMDAGetElementOwnershipRanges3d(DM da,PetscInt **_lx,PetscInt **_ly,PetscInt **_lz)
265: {
266: PetscMPIInt nproc,rank;
267: MPI_Comm comm;
268: PetscInt M,N,P,pM,pN,pP;
269: PetscInt i,j,k,dim,esi,esj,esk,mx,my,mz;
270: PetscInt *lmx,*lmy,*lmz,*tmp;
274: /* create file name */
275: PetscObjectGetComm( (PetscObject)da, &comm );
276: MPI_Comm_size( comm, &nproc );
277: MPI_Comm_rank( comm, &rank );
279: DMDAGetInfo( da, &dim, &M,&N,&P, &pM,&pN,&pP, 0, 0, 0,0,0,0 );
280: DMDAGetElementCorners(da,&esi,&esj,&esk,&mx,&my,&mz);
282: if (dim == 1) {
283: pN = 1;
284: pP = 1;
285: }
286: if (dim == 2) {
287: pP = 1;
288: }
290: PetscMalloc( sizeof(PetscInt)*(nproc), &tmp );
292: PetscMalloc( sizeof(PetscInt)*(pM+1), &lmx );
293: PetscMalloc( sizeof(PetscInt)*(pN+1), &lmy );
294: PetscMalloc( sizeof(PetscInt)*(pP+1), &lmz );
296: if (dim >= 1) {
297: MPI_Allgather ( &mx, 1, MPIU_INT, tmp, 1, MPIU_INT, comm );
298: j = k = 0;
299: for( i=0; i<pM; i++ ) {
300: PetscInt procid = i + j*pM; /* convert proc(i,j,k) to pid */
301: lmx[i] = tmp[procid];
302: }
303: }
305: if (dim >= 2 ) {
306: MPI_Allgather ( &my, 1, MPIU_INT, tmp, 1, MPIU_INT, comm );
307: i = k = 0;
308: for( j=0; j<pN; j++ ) {
309: PetscInt procid = i + j*pM; /* convert proc(i,j,k) to pid */
310: lmy[j] = tmp[procid];
311: }
312: }
314: if (dim == 3 ) {
315: MPI_Allgather ( &mz, 1, MPIU_INT, tmp, 1, MPIU_INT, comm );
316: i = j = 0;
317: for( k=0; k<pP; k++ ) {
318: PetscInt procid = i + j*pM + k*pM*pN; /* convert proc(i,j,k) to pid */
319: lmz[k] = tmp[procid];
320: }
321: }
323: if(_lx) { *_lx = lmx; }
324: if(_ly) { *_ly = lmy; }
325: if(_lz) { *_lz = lmz; }
327: PetscFree(tmp);
329: return(0);
330: }
332: /*
333: i,j are the element indices
334: The unknown is a vector quantity.
335: The s[].c is used to indicate the degree of freedom.
336: */
339: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j,PetscInt k)
340: {
341: PetscInt n;
343: /* velocity */
344: n = 0;
345: /* node 0 */
346: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++; /* Vx0 */
347: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++; /* Vy0 */
348: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++; /* Vz0 */
350: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
351: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
352: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
354: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
355: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
356: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
358: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++;
359: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++;
360: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++;
362: /* */
363: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++; /* Vx4 */
364: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++; /* Vy4 */
365: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++; /* Vz4 */
367: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
368: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
369: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
371: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
372: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
373: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
375: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++;
376: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++;
377: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++;
379: /* pressure */
380: n = 0;
381: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++; /* P0 */
382: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
383: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
384: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++;
386: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++; /* P0 */
387: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
388: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
389: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++;
391: return(0);
392: }
396: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords,PetscInt i,PetscInt j,PetscInt k,PetscScalar el_coord[])
397: {
399: /* get coords for the element */
400: el_coord[0] = coords[k ][j ][i ].x;
401: el_coord[1] = coords[k ][j ][i ].y;
402: el_coord[2] = coords[k ][j ][i ].z;
404: el_coord[3] = coords[k ][j+1][i].x;
405: el_coord[4] = coords[k ][j+1][i].y;
406: el_coord[5] = coords[k ][j+1][i].z;
408: el_coord[6] = coords[k ][j+1][i+1].x;
409: el_coord[7] = coords[k ][j+1][i+1].y;
410: el_coord[8] = coords[k ][j+1][i+1].z;
412: el_coord[9 ] = coords[k ][j ][i+1].x;
413: el_coord[10] = coords[k ][j ][i+1].y;
414: el_coord[11] = coords[k ][j ][i+1].z;
416: el_coord[12] = coords[k+1][j ][i ].x;
417: el_coord[13] = coords[k+1][j ][i ].y;
418: el_coord[14] = coords[k+1][j ][i ].z;
420: el_coord[15] = coords[k+1][j+1][i ].x;
421: el_coord[16] = coords[k+1][j+1][i ].y;
422: el_coord[17] = coords[k+1][j+1][i ].z;
424: el_coord[18] = coords[k+1][j+1][i+1].x;
425: el_coord[19] = coords[k+1][j+1][i+1].y;
426: el_coord[20] = coords[k+1][j+1][i+1].z;
428: el_coord[21] = coords[k+1][j ][i+1].x;
429: el_coord[22] = coords[k+1][j ][i+1].y;
430: el_coord[23] = coords[k+1][j ][i+1].z;
431: return(0);
432: }
436: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field,PetscInt i,PetscInt j,PetscInt k,StokesDOF nodal_fields[])
437: {
439: /* get the nodal fields for u */
440: nodal_fields[0].u_dof = field[k ][j ][i ].u_dof;
441: nodal_fields[0].v_dof = field[k ][j ][i ].v_dof;
442: nodal_fields[0].w_dof = field[k ][j ][i ].w_dof;
444: nodal_fields[1].u_dof = field[k ][j+1][i ].u_dof;
445: nodal_fields[1].v_dof = field[k ][j+1][i ].v_dof;
446: nodal_fields[1].w_dof = field[k ][j+1][i ].w_dof;
448: nodal_fields[2].u_dof = field[k ][j+1][i+1].u_dof;
449: nodal_fields[2].v_dof = field[k ][j+1][i+1].v_dof;
450: nodal_fields[2].w_dof = field[k ][j+1][i+1].w_dof;
452: nodal_fields[3].u_dof = field[k ][j ][i+1].u_dof;
453: nodal_fields[3].v_dof = field[k ][j ][i+1].v_dof;
454: nodal_fields[3].w_dof = field[k ][j ][i+1].w_dof;
456: nodal_fields[4].u_dof = field[k+1][j ][i ].u_dof;
457: nodal_fields[4].v_dof = field[k+1][j ][i ].v_dof;
458: nodal_fields[4].w_dof = field[k+1][j ][i ].w_dof;
460: nodal_fields[5].u_dof = field[k+1][j+1][i ].u_dof;
461: nodal_fields[5].v_dof = field[k+1][j+1][i ].v_dof;
462: nodal_fields[5].w_dof = field[k+1][j+1][i ].w_dof;
464: nodal_fields[6].u_dof = field[k+1][j+1][i+1].u_dof;
465: nodal_fields[6].v_dof = field[k+1][j+1][i+1].v_dof;
466: nodal_fields[6].w_dof = field[k+1][j+1][i+1].w_dof;
468: nodal_fields[7].u_dof = field[k+1][j ][i+1].u_dof;
469: nodal_fields[7].v_dof = field[k+1][j ][i+1].v_dof;
470: nodal_fields[7].w_dof = field[k+1][j ][i+1].w_dof;
472: /* get the nodal fields for p */
473: nodal_fields[0].p_dof = field[k ][j ][i ].p_dof;
474: nodal_fields[1].p_dof = field[k ][j+1][i ].p_dof;
475: nodal_fields[2].p_dof = field[k ][j+1][i+1].p_dof;
476: nodal_fields[3].p_dof = field[k ][j ][i+1].p_dof;
478: nodal_fields[4].p_dof = field[k+1][j ][i ].p_dof;
479: nodal_fields[5].p_dof = field[k+1][j+1][i ].p_dof;
480: nodal_fields[6].p_dof = field[k+1][j+1][i+1].p_dof;
481: nodal_fields[7].p_dof = field[k+1][j ][i+1].p_dof;
483: return(0);
484: }
486: static PetscInt ASS_MAP_wIwDI_uJuDJ(
487: PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,
488: PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)
489: {
490: PetscInt ij;
491: PetscInt r,c,nc;
493: nc = u_NPE*u_dof;
495: r = w_dof*wi+wd;
496: c = u_dof*ui+ud;
498: ij = r*nc+c;
500: return ij;
501: }
505: static PetscErrorCode DMDASetValuesLocalStencil3D_ADD_VALUES(StokesDOF ***fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
506: {
507: PetscInt n,I,J,K;
510: for (n = 0; n<NODES_PER_EL; n++) {
512: I = u_eqn[NSD*n ].i;
513: J = u_eqn[NSD*n ].j;
514: K = u_eqn[NSD*n ].k;
515: fields_F[K][J][I].u_dof = fields_F[K][J][I].u_dof+Fe_u[NSD*n ];
517: I = u_eqn[NSD*n+1].i;
518: J = u_eqn[NSD*n+1].j;
519: K = u_eqn[NSD*n+1].k;
520: fields_F[K][J][I].v_dof = fields_F[K][J][I].v_dof+Fe_u[NSD*n+1];
522: I = u_eqn[NSD*n+2].i;
523: J = u_eqn[NSD*n+2].j;
524: K = u_eqn[NSD*n+2].k;
525: fields_F[K][J][I].w_dof = fields_F[K][J][I].w_dof+Fe_u[NSD*n+2];
528: I = p_eqn[n].i;
529: J = p_eqn[n].j;
530: K = p_eqn[n].k;
531: fields_F[K][J][I].p_dof = fields_F[K][J][I].p_dof+Fe_p[n ];
533: }
534: return(0);
535: }
537: static void FormStressOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
538: {
539: PetscInt ngp;
540: PetscScalar gp_xi[GAUSS_POINTS][NSD];
541: PetscScalar gp_weight[GAUSS_POINTS];
542: PetscInt p,i,j,k;
543: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
544: PetscScalar J_p,tildeD[6];
545: PetscScalar B[6][U_DOFS*NODES_PER_EL];
546: const PetscInt nvdof = U_DOFS*NODES_PER_EL;
548: /* define quadrature rule */
549: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
551: /* evaluate integral */
552: for (p = 0; p < ngp; p++) {
553: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
554: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
556: for (i = 0; i < NODES_PER_EL; i++) {
557: PetscScalar d_dx_i = GNx_p[0][i];
558: PetscScalar d_dy_i = GNx_p[1][i];
559: PetscScalar d_dz_i = GNx_p[2][i];
561: B[0][3*i ] = d_dx_i; B[0][3*i+1] = 0.0; B[0][3*i+2] = 0.0;
562: B[1][3*i ] = 0.0; B[1][3*i+1] = d_dy_i; B[1][3*i+2] = 0.0;
563: B[2][3*i ] = 0.0; B[2][3*i+1] = 0.0; B[2][3*i+2] = d_dz_i;
565: B[3][3*i] = d_dy_i; B[3][3*i+1] = d_dx_i; B[3][3*i+2] = 0.0; /* e_xy */
566: B[4][3*i] = d_dz_i; B[4][3*i+1] = 0.0; B[4][3*i+2] = d_dx_i;/* e_xz */
567: B[5][3*i] = 0.0; B[5][3*i+1] = d_dz_i; B[5][3*i+2] = d_dy_i;/* e_yz */
568: }
571: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
572: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
573: tildeD[2] = 2.0*gp_weight[p]*J_p*eta[p];
575: tildeD[3] = gp_weight[p]*J_p*eta[p];
576: tildeD[4] = gp_weight[p]*J_p*eta[p];
577: tildeD[5] = gp_weight[p]*J_p*eta[p];
579: /* form Bt tildeD B */
580: /*
581: Ke_ij = Bt_ik . D_kl . B_lj
582: = B_ki . D_kl . B_lj
583: = B_ki . D_kk . B_kj
584: */
585: for (i = 0; i < nvdof; i++) {
586: for (j = i; j < nvdof; j++) {
587: for (k = 0; k < 6; k++) {
588: Ke[i*nvdof+j] += B[k][i]*tildeD[k]*B[k][j];
589: }
590: }
591: }
593: }
594: /* fill lower triangular part */
595: #ifdef ASSEMBLE_LOWER_TRIANGULAR
596: for (i = 0; i < nvdof; i++) {
597: for (j = i; j < nvdof; j++) {
598: Ke[j*nvdof+i] = Ke[i*nvdof+j];
599: }
600: }
601: #endif
602: }
604: static void FormGradientOperatorQ13D(PetscScalar Ke[],PetscScalar coords[])
605: {
606: PetscInt ngp;
607: PetscScalar gp_xi[GAUSS_POINTS][NSD];
608: PetscScalar gp_weight[GAUSS_POINTS];
609: PetscInt p,i,j,di;
610: PetscScalar Ni_p[NODES_PER_EL];
611: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
612: PetscScalar J_p,fac;
614: /* define quadrature rule */
615: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
617: /* evaluate integral */
618: for (p = 0; p < ngp; p++) {
619: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
620: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
621: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
622: fac = gp_weight[p]*J_p;
624: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
625: for (di = 0; di < NSD; di++) { /* u dofs */
626: for (j = 0; j < NODES_PER_EL; j++) { /* p nodes, p dofs = 1 (ie no loop) */
627: PetscInt IJ;
628: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,3,j,0,NODES_PER_EL,1);
630: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
631: }
632: }
633: }
634: }
635: }
637: static void FormDivergenceOperatorQ13D(PetscScalar De[],PetscScalar coords[])
638: {
639: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
640: PetscInt i,j;
641: PetscInt nr_g,nc_g;
643: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
644: FormGradientOperatorQ13D(Ge,coords);
646: nr_g = U_DOFS*NODES_PER_EL;
647: nc_g = P_DOFS*NODES_PER_EL;
649: for (i = 0; i < nr_g; i++) {
650: for (j = 0; j < nc_g; j++) {
651: De[nr_g*j+i] = Ge[nc_g*i+j];
652: }
653: }
654: }
656: static void FormStabilisationOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
657: {
658: PetscInt ngp;
659: PetscScalar gp_xi[GAUSS_POINTS][NSD];
660: PetscScalar gp_weight[GAUSS_POINTS];
661: PetscInt p,i,j;
662: PetscScalar Ni_p[NODES_PER_EL];
663: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
664: PetscScalar J_p,fac,eta_avg;
666: /* define quadrature rule */
667: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
669: /* evaluate integral */
670: for (p = 0; p < ngp; p++) {
671: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
672: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
673: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
674: fac = gp_weight[p]*J_p;
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*(Ni_p[i]*Ni_p[j]-0.015625);
679: }
680: }
681: */
683: for (i = 0; i < NODES_PER_EL; i++) {
684: for (j = 0; j < NODES_PER_EL; j++) {
685: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
686: }
687: }
689: }
691: /* scale */
692: eta_avg = 0.0;
693: for (p = 0; p < ngp; p++) {
694: eta_avg += eta[p];
695: }
696: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
697: fac = 1.0/eta_avg;
699: /*
700: for (i = 0; i < NODES_PER_EL; i++) {
701: for (j = i; j < NODES_PER_EL; j++) {
702: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
703: #ifdef ASSEMBLE_LOWER_TRIANGULAR
704: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
705: #endif
706: }
707: }
708: */
710: for (i = 0; i < NODES_PER_EL; i++) {
711: for (j = 0; j < NODES_PER_EL; j++) {
712: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
713: }
714: }
715: }
717: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
718: {
719: PetscInt ngp;
720: PetscScalar gp_xi[GAUSS_POINTS][NSD];
721: PetscScalar gp_weight[GAUSS_POINTS];
722: PetscInt p,i,j;
723: PetscScalar Ni_p[NODES_PER_EL];
724: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
725: PetscScalar J_p,fac,eta_avg;
727: /* define quadrature rule */
728: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
730: /* evaluate integral */
731: for (p = 0; p < ngp; p++) {
732: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
733: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
734: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
735: fac = gp_weight[p]*J_p;
737: /*
738: for (i = 0; i < NODES_PER_EL; i++) {
739: for (j = i; j < NODES_PER_EL; j++) {
740: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
741: }
742: }
743: */
745: for (i = 0; i < NODES_PER_EL; i++) {
746: for (j = 0; j < NODES_PER_EL; j++) {
747: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
748: }
749: }
751: }
753: /* scale */
754: eta_avg = 0.0;
755: for (p = 0; p < ngp; p++) {
756: eta_avg += eta[p];
757: }
758: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
759: fac = 1.0/eta_avg;
760: /*
761: for (i = 0; i < NODES_PER_EL; i++) {
762: for (j = i; j < NODES_PER_EL; j++) {
763: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
764: #ifdef ASSEMBLE_LOWER_TRIANGULAR
765: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
766: #endif
767: }
768: }
769: */
771: for (i = 0; i < NODES_PER_EL; i++) {
772: for (j = 0; j < NODES_PER_EL; j++) {
773: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
774: }
775: }
777: }
779: static void FormMomentumRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[],PetscScalar fz[])
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++) {
800: Fe[NSD*i ] -= fac*Ni_p[i]*fx[p];
801: Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
802: Fe[NSD*i+2] -= fac*Ni_p[i]*fz[p];
803: }
804: }
805: }
807: static void FormContinuityRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar hc[])
808: {
809: PetscInt ngp;
810: PetscScalar gp_xi[GAUSS_POINTS][NSD];
811: PetscScalar gp_weight[GAUSS_POINTS];
812: PetscInt p,i;
813: PetscScalar Ni_p[NODES_PER_EL];
814: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
815: PetscScalar J_p,fac;
817: /* define quadrature rule */
818: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
820: /* evaluate integral */
821: for (p = 0; p < ngp; p++) {
822: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
823: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
824: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
825: fac = gp_weight[p]*J_p;
827: for (i = 0; i < NODES_PER_EL; i++) {
828: Fe[i] -= fac*Ni_p[i]*hc[p];
829: }
830: }
831: }
833: #define _ZERO_ROWCOL_i(A,i) { \
834: PetscInt KK; \
835: PetscScalar tmp = A[24*(i)+(i)]; \
836: for(KK=0;KK<24;KK++){A[24*(i)+KK]=0.0;} \
837: for(KK=0;KK<24;KK++){A[24*KK+(i)]=0.0;} \
838: A[24*(i)+(i)] = tmp;} \
840: #define _ZERO_ROW_i(A,i) { \
841: PetscInt KK; \
842: for(KK=0;KK<8;KK++){A[8*(i)+KK]=0.0;}}
844: #define _ZERO_COL_i(A,i) { \
845: PetscInt KK; \
846: for(KK=0;KK<8;KK++){A[24*KK+(i)]=0.0;}}
850: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
851: {
852: DM cda;
853: Vec coords;
854: DMDACoor3d ***_coords;
855: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
856: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
857: PetscInt sex,sey,sez,mx,my,mz;
858: PetscInt ei,ej,ek;
859: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
860: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
861: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
862: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
863: PetscScalar el_coords[NODES_PER_EL*NSD];
864: Vec local_properties;
865: GaussPointCoefficients ***props;
866: PetscScalar *prop_eta;
867: PetscInt n,M,N,P;
868: PetscLogDouble t0,t1;
869: PetscErrorCode ierr;
873: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
874: /* setup for coords */
875: DMDAGetCoordinateDA(stokes_da,&cda);
876: DMDAGetGhostedCoordinates(stokes_da,&coords);
877: DMDAVecGetArray(cda,coords,&_coords);
879: /* setup for coefficients */
880: DMCreateLocalVector(properties_da,&local_properties);
881: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
882: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
883: DMDAVecGetArray(properties_da,local_properties,&props);
885: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
886: PetscGetTime(&t0);
887: for (ek = sez; ek < sez+mz; ek++) {
888: for (ej = sey; ej < sey+my; ej++) {
889: for (ei = sex; ei < sex+mx; ei++) {
890: /* get coords for the element */
891: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
893: /* get coefficients for the element */
894: prop_eta = props[ek][ej][ei].eta;
896: /* initialise element stiffness matrix */
897: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
898: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
899: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
900: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
902: /* form element stiffness matrix */
903: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
904: FormGradientOperatorQ13D(Ge,el_coords);
905: /*#ifdef ASSEMBLE_LOWER_TRIANGULAR*/
906: FormDivergenceOperatorQ13D(De,el_coords);
907: /*#endif*/
908: FormStabilisationOperatorQ13D(Ce,el_coords,prop_eta);
910: /* insert element matrix into global matrix */
911: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
913: for (n=0; n<NODES_PER_EL; n++) {
914: if ( (u_eqn[3*n ].i == 0) || (u_eqn[3*n ].i == M-1) ) {
915: _ZERO_ROWCOL_i(Ae,3*n);
916: _ZERO_ROW_i (Ge,3*n);
917: _ZERO_COL_i (De,3*n);
918: }
920: if ( (u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1) ) {
921: _ZERO_ROWCOL_i(Ae,3*n+1);
922: _ZERO_ROW_i (Ge,3*n+1);
923: _ZERO_COL_i (De,3*n+1);
924: }
926: if ( (u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1) ) {
927: _ZERO_ROWCOL_i(Ae,3*n+2);
928: _ZERO_ROW_i (Ge,3*n+2);
929: _ZERO_COL_i (De,3*n+2);
930: }
931: }
932: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
933: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
934: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
935: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
936: }
937: }
938: }
939: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
940: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
942: DMDAVecRestoreArray(cda,coords,&_coords);
944: DMDAVecRestoreArray(properties_da,local_properties,&props);
945: VecDestroy(&local_properties);
946: PetscGetTime(&t1);
948: return(0);
949: }
953: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
954: {
955: DM cda;
956: Vec coords;
957: DMDACoor3d ***_coords;
958: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
959: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
960: PetscInt sex,sey,sez,mx,my,mz;
961: PetscInt ei,ej,ek;
962: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
963: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
964: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
965: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
966: PetscScalar el_coords[NODES_PER_EL*NSD];
967: Vec local_properties;
968: GaussPointCoefficients ***props;
969: PetscScalar *prop_eta;
970: PetscInt n,M,N,P;
971: PetscErrorCode ierr;
975: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
976: /* setup for coords */
977: DMDAGetCoordinateDA(stokes_da,&cda);
978: DMDAGetGhostedCoordinates(stokes_da,&coords);
979: DMDAVecGetArray(cda,coords,&_coords);
981: /* setup for coefficients */
982: DMCreateLocalVector(properties_da,&local_properties);
983: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
984: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
985: DMDAVecGetArray(properties_da,local_properties,&props);
987: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
988: for (ek = sez; ek < sez+mz; ek++) {
989: for (ej = sey; ej < sey+my; ej++) {
990: for (ei = sex; ei < sex+mx; ei++) {
991: /* get coords for the element */
992: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
994: /* get coefficients for the element */
995: prop_eta = props[ek][ej][ei].eta;
997: /* initialise element stiffness matrix */
998: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
999: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
1000: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
1001: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
1003: /* form element stiffness matrix */
1004: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
1005: FormGradientOperatorQ13D(Ge,el_coords);
1006: /*FormDivergenceOperatorQ13D(De,el_coords);*/
1007: FormScaledMassMatrixOperatorQ13D(Ce,el_coords,prop_eta);
1009: /* insert element matrix into global matrix */
1010: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
1012: for (n=0; n<NODES_PER_EL; n++) {
1013: if ( (u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1) ) {
1014: _ZERO_ROWCOL_i(Ae,3*n);
1015: _ZERO_ROW_i (Ge,3*n);
1016: _ZERO_COL_i (De,3*n);
1017: }
1019: if ( (u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1) ) {
1020: _ZERO_ROWCOL_i(Ae,3*n+1);
1021: _ZERO_ROW_i (Ge,3*n+1);
1022: _ZERO_COL_i (De,3*n+1);
1023: }
1025: if ( (u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1) ) {
1026: _ZERO_ROWCOL_i(Ae,3*n+2);
1027: _ZERO_ROW_i (Ge,3*n+2);
1028: _ZERO_COL_i (De,3*n+2);
1029: }
1030: }
1031: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
1032: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
1033: /*MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);*/
1034: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
1035: }
1036: }
1037: }
1038: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1039: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1041: DMDAVecRestoreArray(cda,coords,&_coords);
1043: DMDAVecRestoreArray(properties_da,local_properties,&props);
1044: VecDestroy(&local_properties);
1045: return(0);
1046: }
1050: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
1051: {
1052: DM cda;
1053: Vec coords;
1054: DMDACoor3d ***_coords;
1055: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
1056: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
1057: PetscInt sex,sey,sez,mx,my,mz;
1058: PetscInt ei,ej,ek;
1059: PetscScalar Fe[NODES_PER_EL*U_DOFS];
1060: PetscScalar He[NODES_PER_EL*P_DOFS];
1061: PetscScalar el_coords[NODES_PER_EL*NSD];
1062: Vec local_properties;
1063: GaussPointCoefficients ***props;
1064: PetscScalar *prop_fx,*prop_fy,*prop_fz,*prop_hc;
1065: Vec local_F;
1066: StokesDOF ***ff;
1067: PetscInt n,M,N,P;
1068: PetscErrorCode ierr;
1072: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
1073: /* setup for coords */
1074: DMDAGetCoordinateDA(stokes_da,&cda);
1075: DMDAGetGhostedCoordinates(stokes_da,&coords);
1076: DMDAVecGetArray(cda,coords,&_coords);
1078: /* setup for coefficients */
1079: DMGetLocalVector(properties_da,&local_properties);
1080: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
1081: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
1082: DMDAVecGetArray(properties_da,local_properties,&props);
1084: /* get acces to the vector */
1085: DMGetLocalVector(stokes_da,&local_F);
1086: VecZeroEntries(local_F);
1087: DMDAVecGetArray(stokes_da,local_F,&ff);
1090: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
1091: for (ek = sez; ek < sez+mz; ek++) {
1092: for (ej = sey; ej < sey+my; ej++) {
1093: for (ei = sex; ei < sex+mx; ei++) {
1094: /* get coords for the element */
1095: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1097: /* get coefficients for the element */
1098: prop_fx = props[ek][ej][ei].fx;
1099: prop_fy = props[ek][ej][ei].fy;
1100: prop_fz = props[ek][ej][ei].fz;
1101: prop_hc = props[ek][ej][ei].hc;
1103: /* initialise element stiffness matrix */
1104: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
1105: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
1107: /* form element stiffness matrix */
1108: FormMomentumRhsQ13D(Fe,el_coords,prop_fx,prop_fy,prop_fz);
1109: FormContinuityRhsQ13D(He,el_coords,prop_hc);
1111: /* insert element matrix into global matrix */
1112: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
1115: for (n=0; n<NODES_PER_EL; n++) {
1116: if ( (u_eqn[3*n ].i == 0) || (u_eqn[3*n ].i == M-1) ) {
1117: Fe[3*n ] = 0.0;
1118: }
1120: if ( (u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1) ) {
1121: Fe[3*n+1] = 0.0;
1122: }
1124: if ( (u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1) ) {
1125: Fe[3*n+2] = 0.0;
1126: }
1127: }
1129: DMDASetValuesLocalStencil3D_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
1130: }
1131: }
1132: }
1133: DMDAVecRestoreArray(stokes_da,local_F,&ff);
1134: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
1135: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
1136: DMRestoreLocalVector(stokes_da,&local_F);
1139: DMDAVecRestoreArray(cda,coords,&_coords);
1141: DMDAVecRestoreArray(properties_da,local_properties,&props);
1142: DMRestoreLocalVector(properties_da,&local_properties);
1143: return(0);
1144: }
1146: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta,PetscReal *MX,PetscReal *MY,PetscReal *MZ)
1147: {
1148: *theta = 0.0;
1149: *MX = 2.0 * M_PI;
1150: *MY = 2.0 * M_PI;
1151: *MZ = 2.0 * M_PI;
1152: }
1153: static void evaluate_MS_FrankKamentski(PetscReal pos[],PetscReal v[],PetscReal *p,PetscReal *eta,PetscReal Fm[],PetscReal *Fc)
1154: {
1155: PetscReal x,y,z;
1156: PetscReal theta,MX,MY,MZ;
1158: evaluate_MS_FrankKamentski_constants(&theta,&MX,&MY,&MZ);
1159: x = pos[0];
1160: y = pos[1];
1161: z = pos[2];
1162: if (v) {
1163: /*
1164: v[0] = pow(z,3)*exp(y)*sin(M_PI*x);
1165: v[1] = z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y);
1166: v[2] = -(pow(x,3) + pow(y,3))*sin(2.0*M_PI*z);
1167: */
1168: /*
1169: v[0] = sin(M_PI*x);
1170: v[1] = sin(M_PI*y);
1171: v[2] = sin(M_PI*z);
1172: */
1173: v[0] = pow(z,3)*exp(y)*sin(M_PI*x);
1174: v[1] = z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y);
1175: v[2] = pow(z,2)*(cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 - M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2) - M_PI*pow(z,4)*cos(M_PI*x)*exp(y)/4 ;
1176: }
1177: if (p) {
1178: *p = pow(x,2) + pow(y,2) + pow(z,2);
1179: }
1180: if (eta) {
1181: /**eta = exp(-theta*(1.0 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)));*/
1182: *eta = 1.0;
1183: }
1184: if (Fm) {
1185: /*
1186: 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))) ;
1187: 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)));
1188: 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))) ;
1189: */
1190: /*
1191: Fm[0]=-2*x - 2.0*pow(M_PI,2)*sin(M_PI*x);
1192: Fm[1]=-2*y - 2.0*pow(M_PI,2)*sin(M_PI*y);
1193: Fm[2]=-2*z - 2.0*pow(M_PI,2)*sin(M_PI*z);
1194: */
1195: /*
1196: 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) ;
1197: 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);
1198: 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) ;
1199: */
1200: Fm[0] = -2*x + 2*z*(pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x) - 1.0*M_PI*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x)) + pow(z,3)*exp(y)*sin(M_PI*x) + 6.0*z*exp(y)*sin(M_PI*x) - 1.0*pow(M_PI,2)*pow(z,3)*exp(y)*sin(M_PI*x) + 2.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x) - 2.0*z*pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x);
1201: Fm[1] = -2*y + 2*z*(-cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)) + 2.0*z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 6.0*z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 4.0*M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;
1202: Fm[2] = -2*z + pow(z,2)*(-2.0*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 2.0*pow(M_PI,3)*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)) + pow(z,2)*(cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 - 3*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + pow(M_PI,3)*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2 - 3*M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2) + 1.0*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 0.25*pow(M_PI,3)*pow(z,4)*cos(M_PI*x)*exp(y) - 0.25*M_PI*pow(z,4)*cos(M_PI*x)*exp(y) - 3.0*M_PI*pow(z,2)*cos(M_PI*x)*exp(y) - 1.0*M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;
1203: }
1204: if (Fc) {
1205: /**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) ;*/
1206: /**Fc = M_PI*cos(M_PI*x) + M_PI*cos(M_PI*y) + M_PI*cos(M_PI*z);*/
1207: /**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);*/
1208: *Fc = 0.0;
1209: }
1210: }
1214: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx,PetscInt my,PetscInt mz,DM *_da,Vec *_X)
1215: {
1216: DM da,cda;
1217: Vec X,local_X;
1218: StokesDOF ***_stokes;
1219: Vec coords;
1220: DMDACoor3d ***_coords;
1221: PetscInt si,sj,sk,ei,ej,ek,i,j,k;
1225: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1226: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,4,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&da);
1227: DMDASetFieldName(da,0,"anlytic_Vx");
1228: DMDASetFieldName(da,1,"anlytic_Vy");
1229: DMDASetFieldName(da,2,"anlytic_Vz");
1230: DMDASetFieldName(da,3,"analytic_P");
1232: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
1234: DMDAGetGhostedCoordinates(da,&coords);
1235: DMDAGetCoordinateDA(da,&cda);
1236: DMDAVecGetArray(cda,coords,&_coords);
1238: DMCreateGlobalVector(da,&X);
1239: DMCreateLocalVector(da,&local_X);
1240: DMDAVecGetArray(da,local_X,&_stokes);
1242: DMDAGetGhostCorners(da,&si,&sj,&sk,&ei,&ej,&ek);
1243: for (k = sk; k < sk+ek; k++) {
1244: for (j = sj; j < sj+ej; j++) {
1245: for (i = si; i < si+ei; i++) {
1246: PetscReal pos[NSD],pressure,vel[NSD];
1248: pos[0] = PetscRealPart(_coords[k][j][i].x);
1249: pos[1] = PetscRealPart(_coords[k][j][i].y);
1250: pos[2] = PetscRealPart(_coords[k][j][i].z);
1252: evaluate_MS_FrankKamentski(pos,vel,&pressure,PETSC_NULL,PETSC_NULL,PETSC_NULL);
1254: _stokes[k][j][i].u_dof = vel[0];
1255: _stokes[k][j][i].v_dof = vel[1];
1256: _stokes[k][j][i].w_dof = vel[2];
1257: _stokes[k][j][i].p_dof = pressure;
1258: }
1259: }
1260: }
1261: DMDAVecRestoreArray(da,local_X,&_stokes);
1262: DMDAVecRestoreArray(cda,coords,&_coords);
1264: DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1265: DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);
1267: VecDestroy(&local_X);
1269: *_da = da;
1270: *_X = X;
1272: return(0);
1273: }
1277: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da,Vec X,Vec X_analytic)
1278: {
1279: DM cda;
1280: Vec coords,X_analytic_local,X_local;
1281: DMDACoor3d ***_coords;
1282: PetscInt sex,sey,sez,mx,my,mz;
1283: PetscInt ei,ej,ek;
1284: PetscScalar el_coords[NODES_PER_EL*NSD];
1285: StokesDOF ***stokes_analytic,***stokes;
1286: StokesDOF stokes_analytic_e[NODES_PER_EL],stokes_e[NODES_PER_EL];
1288: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1289: PetscScalar Ni_p[NODES_PER_EL];
1290: PetscInt ngp;
1291: PetscScalar gp_xi[GAUSS_POINTS][NSD];
1292: PetscScalar gp_weight[GAUSS_POINTS];
1293: PetscInt p,i;
1294: PetscScalar J_p,fac;
1295: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1296: PetscScalar tint_p_ms,tint_p,int_p_ms,int_p;
1297: PetscInt M;
1298: PetscReal xymin[NSD],xymax[NSD];
1302: /* define quadrature rule */
1303: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1305: /* setup for coords */
1306: DMDAGetCoordinateDA(stokes_da,&cda);
1307: DMDAGetGhostedCoordinates(stokes_da,&coords);
1308: DMDAVecGetArray(cda,coords,&_coords);
1310: /* setup for analytic */
1311: DMCreateLocalVector(stokes_da,&X_analytic_local);
1312: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1313: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1314: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1316: /* setup for solution */
1317: DMCreateLocalVector(stokes_da,&X_local);
1318: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1319: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1320: DMDAVecGetArray(stokes_da,X_local,&stokes);
1322: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1323: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1325: h = (xymax[0]-xymin[0])/((PetscReal)(M-1));
1327: tp_L2 = tu_L2 = tu_H1 = 0.0;
1328: tint_p_ms = tint_p = 0.0;
1330: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
1332: for (ek = sez; ek < sez+mz; ek++) {
1333: for (ej = sey; ej < sey+my; ej++) {
1334: for (ei = sex; ei < sex+mx; ei++) {
1335: /* get coords for the element */
1336: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1337: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1338: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1340: /* evaluate integral */
1341: for (p = 0; p < ngp; p++) {
1342: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1343: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1344: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1345: fac = gp_weight[p]*J_p;
1347: for (i = 0; i < NODES_PER_EL; i++) {
1348: tint_p_ms = tint_p_ms+fac*Ni_p[i]*stokes_analytic_e[i].p_dof;
1349: tint_p = tint_p +fac*Ni_p[i]*stokes_e[i].p_dof;
1350: }
1351: }
1352: }
1353: }
1354: }
1355: MPI_Allreduce(&tint_p_ms,&int_p_ms,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1356: MPI_Allreduce(&tint_p,&int_p,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1358: PetscPrintf(PETSC_COMM_WORLD,"\\int P dv %1.4e (h) %1.4e (ms)\n",PetscRealPart(int_p),PetscRealPart(int_p_ms));
1360: /* remove mine and add manufacture one */
1361: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1362: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1364: {
1365: PetscInt k,L,dof;
1366: PetscScalar *fields;
1367: DMDAGetInfo(stokes_da,0, 0,0,0, 0,0,0, &dof,0,0,0,0,0);
1369: VecGetLocalSize(X_local,&L);
1370: VecGetArray(X_local,&fields);
1371: for (k=0; k<L/dof; k++) {
1372: fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1373: }
1374: VecRestoreArray(X_local,&fields);
1376: VecGetLocalSize(X,&L);
1377: VecGetArray(X,&fields);
1378: for (k=0; k<L/dof; k++) {
1379: fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1380: }
1381: VecRestoreArray(X,&fields);
1382: }
1385: DMDAVecGetArray(stokes_da,X_local,&stokes);
1386: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1388: for (ek = sez; ek < sez+mz; ek++) {
1389: for (ej = sey; ej < sey+my; ej++) {
1390: for (ei = sex; ei < sex+mx; ei++) {
1391: /* get coords for the element */
1392: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1393: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1394: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1396: /* evaluate integral */
1397: p_e_L2 = 0.0;
1398: u_e_L2 = 0.0;
1399: u_e_H1 = 0.0;
1400: for (p = 0; p < ngp; p++) {
1401: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1402: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1403: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1404: fac = gp_weight[p]*J_p;
1406: for (i = 0; i < NODES_PER_EL; i++) {
1407: PetscScalar u_error,v_error,w_error;
1409: 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);
1411: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1412: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1413: w_error = stokes_e[i].w_dof-stokes_analytic_e[i].w_dof;
1414: /*
1415: if (p==0) {
1416: printf("p=0: %d %d %d %1.4e,%1.4e,%1.4e \n", ei,ej,ek,u_error,v_error,w_error);
1417: }
1418: */
1419: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error+w_error*w_error);
1421: u_e_H1 = u_e_H1+fac*( GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1422: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error
1423: +GNx_p[2][i]*u_error*GNx_p[2][i]*u_error
1424: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1425: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error
1426: +GNx_p[2][i]*v_error*GNx_p[2][i]*v_error
1427: +GNx_p[0][i]*w_error*GNx_p[0][i]*w_error /* dw/dx */
1428: +GNx_p[1][i]*w_error*GNx_p[1][i]*w_error
1429: +GNx_p[2][i]*w_error*GNx_p[2][i]*w_error );
1430: }
1431: }
1433: tp_L2 += p_e_L2;
1434: tu_L2 += u_e_L2;
1435: tu_H1 += u_e_H1;
1436: }
1437: }
1438: }
1439: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1440: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1441: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1442: p_L2 = PetscSqrtScalar(p_L2);
1443: u_L2 = PetscSqrtScalar(u_L2);
1444: u_H1 = PetscSqrtScalar(u_H1);
1446: PetscPrintf(PETSC_COMM_WORLD,"%1.4e %1.4e %1.4e %1.4e \n",PetscRealPart(h),PetscRealPart(p_L2),PetscRealPart(u_L2),PetscRealPart(u_H1));
1449: DMDAVecRestoreArray(cda,coords,&_coords);
1451: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1452: VecDestroy(&X_analytic_local);
1453: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1454: VecDestroy(&X_local);
1455: return(0);
1456: }
1460: PetscErrorCode DAView_3DVTK_StructuredGrid_appended( DM da, Vec FIELD, const char file_prefix[] )
1461: {
1462: char vtk_filename[PETSC_MAX_PATH_LEN];
1463: PetscMPIInt rank;
1464: MPI_Comm comm;
1465: FILE *vtk_fp;
1466: PetscInt si,sj,sk, nx,ny,nz, i;
1467: PetscInt f, n_fields, N;
1468: DM cda;
1469: Vec coords;
1470: Vec l_FIELD;
1471: PetscScalar *_L_FIELD;
1472: PetscInt memory_offset;
1473: PetscScalar *buffer;
1474: PetscLogDouble t0,t1;
1478: PetscGetTime(&t0);
1480: /* create file name */
1481: PetscObjectGetComm( (PetscObject)da, &comm );
1482: MPI_Comm_rank( comm, &rank );
1483: PetscSNPrintf(vtk_filename, sizeof vtk_filename, "subdomain-%s-p%1.4d.vts", file_prefix, rank );
1485: /* open file and write header */
1486: vtk_fp = fopen( vtk_filename, "w" );
1487: if( vtk_fp == NULL ) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"ERROR(DAView_3DVTK_StructuredGrid_appended): Cannot open file = %s \n", vtk_filename );
1489: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<?xml version=\"1.0\"?>\n");
1491: /* coords */
1492: DMDAGetGhostCorners( da, &si,&sj,&sk, &nx,&ny,&nz );
1493: N = nx * ny * nz;
1495: #ifdef PETSC_WORDS_BIGENDIAN
1496: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1497: #else
1498: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1499: #endif
1500: 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);
1501: 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);
1503: memory_offset = 0;
1505: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <CellData></CellData>\n");
1507: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <Points>\n");
1509: /* copy coordinates */
1510: DMDAGetCoordinateDA( da, &cda);
1511: DMDAGetGhostedCoordinates( da,&coords );
1512: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%d\" />\n",memory_offset);
1513: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N*3;
1515: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </Points>\n");
1518: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PointData Scalars=\" ");
1519: DMDAGetInfo( da, 0, 0,0,0, 0,0,0, &n_fields, 0,0,0, 0, 0 );
1520: for( f=0; f<n_fields; f++ ) {
1521: const char *field_name;
1522: DMDAGetFieldName( da, f, &field_name );
1523: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "%s ", field_name );
1524: }
1525: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "\">\n");
1527: for( f=0; f<n_fields; f++ ) {
1528: const char *field_name;
1530: DMDAGetFieldName( da, f, &field_name );
1531: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <DataArray type=\"Float64\" Name=\"%s\" format=\"appended\" offset=\"%d\"/>\n", field_name,memory_offset );
1532: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N;
1533: }
1535: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </PointData>\n");
1536: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </Piece>\n");
1537: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </StructuredGrid>\n");
1539: PetscMalloc( sizeof(PetscScalar)*N, &buffer );
1540: DMGetLocalVector( da, &l_FIELD );
1541: DMGlobalToLocalBegin( da, FIELD, INSERT_VALUES, l_FIELD );
1542: DMGlobalToLocalEnd( da, FIELD, INSERT_VALUES, l_FIELD );
1543: VecGetArray( l_FIELD, &_L_FIELD );
1545: PetscFPrintf(PETSC_COMM_SELF, vtk_fp," <AppendedData encoding=\"raw\">\n");
1546: PetscFPrintf(PETSC_COMM_SELF, vtk_fp,"_");
1548: /* write coordinates */
1549: {
1550: int length = sizeof(PetscScalar)*N*3;
1551: PetscScalar *allcoords;
1553: fwrite( &length,sizeof(int),1,vtk_fp);
1554: VecGetArray(coords,&allcoords);
1555: fwrite( allcoords, sizeof(PetscScalar),3*N, vtk_fp );
1556: VecRestoreArray(coords,&allcoords);
1557: }
1558: /* write fields */
1559: for( f=0; f<n_fields; f++ ) {
1560: int length = sizeof(PetscScalar)*N;
1561: fwrite( &length,sizeof(int),1,vtk_fp);
1562: /* load */
1563: for( i=0; i<N; i++ ) {
1564: buffer[i] = _L_FIELD[ n_fields*i + f ];
1565: }
1567: /* write */
1568: fwrite( buffer,sizeof(PetscScalar),N,vtk_fp);
1569: }
1570: PetscFPrintf(PETSC_COMM_SELF, vtk_fp,"\n </AppendedData>\n");
1572: PetscFree(buffer);
1573: VecRestoreArray( l_FIELD, &_L_FIELD );
1574: DMRestoreLocalVector( da, &l_FIELD );
1576: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "</VTKFile>\n");
1579: if( vtk_fp!= NULL ) {
1580: fclose( vtk_fp );
1581: vtk_fp = NULL;
1582: }
1584: PetscGetTime(&t1);
1585: return(0);
1586: }
1590: PetscErrorCode DAViewVTK_write_PieceExtend( FILE *vtk_fp, PetscInt indent_level, DM da, const char local_file_prefix[] )
1591: {
1592: PetscMPIInt nproc,rank;
1593: MPI_Comm comm;
1594: const PetscInt *lx,*ly,*lz;
1595: PetscInt M,N,P,pM,pN,pP,sum,*olx,*oly,*olz;
1596: PetscInt *osx,*osy,*osz,*oex,*oey,*oez;
1597: PetscInt i,j,k,II,stencil;
1601: /* create file name */
1602: PetscObjectGetComm( (PetscObject)da, &comm );
1603: MPI_Comm_size( comm, &nproc );
1604: MPI_Comm_rank( comm, &rank );
1606: DMDAGetInfo( da, 0, &M,&N,&P, &pM,&pN,&pP, 0, &stencil, 0,0,0, 0 );
1607: DMDAGetOwnershipRanges( da,&lx,&ly,&lz );
1609: /* generate start,end list */
1610: PetscMalloc( sizeof(PetscInt)*(pM+1), &olx );
1611: PetscMalloc( sizeof(PetscInt)*(pN+1), &oly );
1612: PetscMalloc( sizeof(PetscInt)*(pP+1), &olz );
1613: sum = 0;
1614: for( i=0; i<pM; i++ ) {
1615: olx[i] = sum;
1616: sum = sum + lx[i];
1617: } olx[pM] = sum;
1618: sum = 0;
1619: for( i=0; i<pN; i++ ) {
1620: oly[i] = sum;
1621: sum = sum + ly[i];
1622: } oly[pN] = sum;
1623: sum = 0;
1624: for( i=0; i<pP; i++ ) {
1625: olz[i] = sum;
1626: sum = sum + lz[i];
1627: } olz[pP] = sum;
1628: PetscMalloc( sizeof(PetscInt)*(pM), &osx );
1629: PetscMalloc( sizeof(PetscInt)*(pN), &osy );
1630: PetscMalloc( sizeof(PetscInt)*(pP), &osz );
1631: PetscMalloc( sizeof(PetscInt)*(pM), &oex );
1632: PetscMalloc( sizeof(PetscInt)*(pN), &oey );
1633: PetscMalloc( sizeof(PetscInt)*(pP), &oez );
1634: for( i=0; i<pM; i++ ) {
1635: osx[i] = olx[i] - stencil;
1636: oex[i] = olx[i] + lx[i] + stencil;
1637: if(osx[i]<0)osx[i]=0;
1638: if(oex[i]>M)oex[i]=M;
1639: }
1641: for( i=0; i<pN; i++ ) {
1642: osy[i] = oly[i] - stencil;
1643: oey[i] = oly[i] + ly[i] + stencil;
1644: if(osy[i]<0)osy[i]=0;
1645: if(oey[i]>M)oey[i]=N;
1646: }
1647: for( i=0; i<pP; i++ ) {
1648: osz[i] = olz[i] - stencil;
1649: oez[i] = olz[i] + lz[i] + stencil;
1650: if(osz[i]<0)osz[i]=0;
1651: if(oez[i]>P)oez[i]=P;
1652: }
1654: for( k=0;k<pP;k++ ) {
1655: for( j=0;j<pN;j++ ) {
1656: for( i=0;i<pM;i++ ) {
1657: char name[PETSC_MAX_PATH_LEN];
1658: PetscInt procid = i + j*pM + k*pM*pN; /* convert proc(i,j,k) to pid */
1659: PetscSNPrintf(name, sizeof name, "subdomain-%s-p%1.4d.vts", local_file_prefix, procid );
1660: for( II=0; II<indent_level; II++ ) {
1661: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," ");
1662: }
1663: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<Piece Extent=\"%d %d %d %d %d %d\" Source=\"%s\"/>\n",
1664: osx[i],oex[i]-1,
1665: osy[j],oey[j]-1,
1666: osz[k],oez[k]-1, name);
1667: }
1668: }
1669: }
1670: PetscFree( olx );
1671: PetscFree( oly );
1672: PetscFree( olz );
1673: PetscFree( osx );
1674: PetscFree( osy );
1675: PetscFree( osz );
1676: PetscFree( oex );
1677: PetscFree( oey );
1678: PetscFree( oez );
1679: return(0);
1680: }
1684: PetscErrorCode DAView_3DVTK_PStructuredGrid( DM da, const char file_prefix[], const char local_file_prefix[] )
1685: {
1686: MPI_Comm comm;
1687: PetscMPIInt nproc,rank;
1688: char vtk_filename[PETSC_MAX_PATH_LEN];
1689: FILE *vtk_fp;
1690: PetscInt M,N,P, si,sj,sk, nx,ny,nz;
1691: PetscInt i,dofs;
1695: /* only master generates this file */
1696: PetscObjectGetComm( (PetscObject)da, &comm );
1697: MPI_Comm_size( comm, &nproc );
1698: MPI_Comm_rank( comm, &rank );
1700: if( rank != 0 ) { return(0); }
1702: /* create file name */
1703: PetscSNPrintf(vtk_filename,sizeof vtk_filename, "%s.pvts", file_prefix );
1704: vtk_fp = fopen( vtk_filename, "w" );
1705: if( vtk_fp == NULL ) {
1706: printf("ERROR(DAView_3DVTK_PStructuredGrid): Cannot open file = %s \n", vtk_filename );
1707: exit(0);
1708: }
1710: /* (VTK) generate pvts header */
1711: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<?xml version=\"1.0\"?>\n");
1713: #ifdef PETSC_WORDS_BIGENDIAN
1714: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1715: #else
1716: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1717: #endif
1719: /* define size of the nodal mesh based on the cell DM */
1720: DMDAGetInfo( da, 0, &M,&N,&P, 0,0,0, &dofs,0,0,0,0,0 );
1721: DMDAGetGhostCorners( da, &si,&sj,&sk, &nx,&ny,&nz );
1722: 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 */
1725: /* DUMP THE CELL REFERENCES */
1726: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PCellData>\n");
1727: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </PCellData>\n");
1729: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PPoints>\n");
1730: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PDataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"%d\"/>\n",NSD);
1731: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </PPoints>\n");
1733: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PPointData>\n");
1734: for(i=0;i<dofs;i++){
1735: const char *fieldname;
1736: DMDAGetFieldName(da,i,&fieldname);
1737: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n",fieldname);
1738: }
1739: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </PPointData>\n");
1741: /* write out the parallel information */
1742: DAViewVTK_write_PieceExtend(vtk_fp,2,da,local_file_prefix);
1745: /* close the file */
1746: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </PStructuredGrid>\n");
1747: PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "</VTKFile>\n");
1749: if(vtk_fp!=NULL){
1750: fclose( vtk_fp );
1751: vtk_fp = NULL;
1752: }
1754: return(0);
1755: }
1759: PetscErrorCode DAView3DPVTS(DM da, Vec x,const char NAME[])
1760: {
1761: char vts_filename[PETSC_MAX_PATH_LEN];
1762: char pvts_filename[PETSC_MAX_PATH_LEN];
1766: PetscSNPrintf(vts_filename,sizeof vts_filename, "%s-mesh", NAME );
1767: DAView_3DVTK_StructuredGrid_appended(da,x,vts_filename);
1769: PetscSNPrintf(pvts_filename,sizeof pvts_filename, "%s-mesh", NAME );
1770: DAView_3DVTK_PStructuredGrid( da, pvts_filename, vts_filename );
1771: return(0);
1772: }
1776: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
1777: {
1779: PetscReal norms[4];
1780: Vec Br, v,w;
1781: Mat A;
1784: KSPGetOperators( ksp, &A, 0, 0 );
1785: MatGetVecs( A, &w,&v );
1786: VecSetBlockSize(v,4);
1788: KSPBuildResidual( ksp, v, w, &Br );
1790: VecStrideNorm(Br,0,NORM_2,&norms[0]);
1791: VecStrideNorm(Br,1,NORM_2,&norms[1]);
1792: VecStrideNorm(Br,2,NORM_2,&norms[2]);
1793: VecStrideNorm(Br,3,NORM_2,&norms[3]);
1795: VecDestroy(&v);
1796: VecDestroy(&w);
1798: PetscPrintf(PETSC_COMM_WORLD," %d 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]);
1800: return(0);
1801: }
1805: static PetscErrorCode PCMGSetupViaCoarsen(PC pc,DM da_fine)
1806: {
1807: PetscInt nlevels,k;
1808: DM *da_list,*daclist;
1809: Mat R;
1813: nlevels = 1;
1814: PetscOptionsGetInt( PETSC_NULL, "-levels", &nlevels, 0 );
1816: PetscMalloc( sizeof(DM)*nlevels, &da_list );
1817: for( k=0; k<nlevels; k++ ) { da_list[k] = PETSC_NULL; }
1818: PetscMalloc( sizeof(DM)*nlevels, &daclist );
1819: for( k=0; k<nlevels; k++ ) { daclist[k] = PETSC_NULL; }
1821: /* finest grid is nlevels - 1 */
1822: daclist[0] = da_fine;
1823: PetscObjectReference( (PetscObject)da_fine );
1824: DMCoarsenHierarchy(da_fine,nlevels-1,&daclist[1]);
1825: for( k=0; k<nlevels; k++ ) {
1826: da_list[k] = daclist[nlevels-1-k];
1827: DMDASetUniformCoordinates(da_list[k], 0.0,1.0, 0.0,1.0, 0.0,1.0);
1828: }
1830: PCMGSetLevels(pc,nlevels,PETSC_NULL);
1831: PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
1832: PCMGSetGalerkin(pc,PETSC_TRUE);
1834: for( k=1; k<nlevels; k++ ){
1835: DMGetInterpolation(da_list[k-1],da_list[k],&R,PETSC_NULL);
1836: PCMGSetInterpolation(pc,k,R);
1837: MatDestroy(&R);
1838: }
1840: /* tidy up */
1841: for( k=0; k<nlevels; k++ ) {
1842: DMDestroy(&da_list[k]);
1843: }
1844: PetscFree(da_list);
1845: PetscFree(daclist);
1847: return(0);
1848: }
1852: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx,PetscInt my,PetscInt mz)
1853: {
1854: DM da_Stokes,da_prop;
1855: PetscInt u_dof,p_dof,dof,stencil_width;
1856: Mat A,B;
1857: PetscInt mxl,myl,mzl;
1858: DM prop_cda,vel_cda;
1859: Vec prop_coords,vel_coords;
1860: PetscInt si,sj,sk,nx,ny,nz,p;
1861: Vec f,X;
1862: PetscInt prop_dof,prop_stencil_width;
1863: Vec properties,l_properties;
1864: PetscReal dx,dy,dz;
1865: PetscInt M,N,P;
1866: DMDACoor3d ***_prop_coords,***_vel_coords;
1867: GaussPointCoefficients ***element_props;
1868: PetscInt its;
1869: KSP ksp_S;
1870: PetscInt model_definition = 0;
1871: PetscInt cpu_x,cpu_y,cpu_z,*lx = PETSC_NULL,*ly = PETSC_NULL,*lz = PETSC_NULL;
1872: PetscInt ei,ej,ek,sex,sey,sez,Imx,Jmy,Kmz;
1873: PetscErrorCode ierr;
1876: /* Generate the da for velocity and pressure */
1877: /* Num nodes in each direction is mx+1, my+1, mz+1 */
1878: u_dof = U_DOFS; /* Vx, Vy - velocities */
1879: p_dof = P_DOFS; /* p - pressure */
1880: dof = u_dof+p_dof;
1881: stencil_width = 1;
1882: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1883: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,PETSC_NULL,PETSC_NULL,PETSC_NULL,&da_Stokes);
1884: DMDASetFieldName(da_Stokes,0,"Vx");
1885: DMDASetFieldName(da_Stokes,1,"Vy");
1886: DMDASetFieldName(da_Stokes,2,"Vz");
1887: DMDASetFieldName(da_Stokes,3,"P");
1889: /* unit box [0,1] x [0,1] x [0,1] */
1890: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.0,1.0);
1892: /* local number of elements */
1893: DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,&mzl);
1895: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1896: DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,&cpu_z,0,0,0,0,0,0);
1897: DMDAGetElementOwnershipRanges3d(da_Stokes,&lx,&ly,&lz);
1899: prop_dof = (PetscInt)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1900: prop_stencil_width = 0;
1901: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1902: mx,my,mz,cpu_x,cpu_y,cpu_z,prop_dof,prop_stencil_width,lx,ly,lz,&da_prop);
1903: PetscFree(lx);
1904: PetscFree(ly);
1905: PetscFree(lz);
1907: /* define centroid positions */
1908: DMDAGetInfo(da_prop,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
1909: dx = 1.0/((PetscReal)(M));
1910: dy = 1.0/((PetscReal)(N));
1911: dz = 1.0/((PetscReal)(P));
1913: DMDASetUniformCoordinates(da_prop,0.0+0.5*dx,1.0-0.5*dx,0.0+0.5*dy,1.0-0.5*dy,0.0+0.5*dz,1.0-0.5*dz);
1915: DMCreateGlobalVector(da_prop,&properties);
1916: DMCreateLocalVector(da_prop,&l_properties);
1917: DMDAVecGetArray(da_prop,l_properties,&element_props);
1919: DMDAGetCoordinateDA(da_prop,&prop_cda);
1920: DMDAGetGhostedCoordinates(da_prop,&prop_coords);
1921: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
1923: DMDAGetGhostCorners(prop_cda,&si,&sj,&sk,&nx,&ny,&nz);
1925: DMDAGetCoordinateDA(da_Stokes,&vel_cda);
1926: DMDAGetGhostedCoordinates(da_Stokes,&vel_coords);
1927: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1930: /* interpolate the coordinates */
1931: DMDAGetElementCorners(da_Stokes,&sex,&sey,&sez,&Imx,&Jmy,&Kmz);
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: /* get coords for the element */
1936: PetscInt ngp;
1937: PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
1938: PetscScalar el_coords[NSD*NODES_PER_EL];
1940: GetElementCoords3D(_vel_coords,ei,ej,ek,el_coords);
1941: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1943: for (p = 0; p < GAUSS_POINTS; p++) {
1944: PetscScalar xi_p[NSD],Ni_p[NODES_PER_EL];
1945: PetscScalar gp_x,gp_y,gp_z;
1946: PetscInt n;
1948: xi_p[0] = gp_xi[p][0];
1949: xi_p[1] = gp_xi[p][1];
1950: xi_p[2] = gp_xi[p][2];
1951: ShapeFunctionQ13D_Evaluate(xi_p,Ni_p);
1953: gp_x = gp_y = gp_z = 0.0;
1954: for (n = 0; n < NODES_PER_EL; n++) {
1955: gp_x = gp_x+Ni_p[n]*el_coords[NSD*n ];
1956: gp_y = gp_y+Ni_p[n]*el_coords[NSD*n+1];
1957: gp_z = gp_z+Ni_p[n]*el_coords[NSD*n+2];
1958: }
1959: element_props[ek][ej][ei].gp_coords[NSD*p ] = gp_x;
1960: element_props[ek][ej][ei].gp_coords[NSD*p+1] = gp_y;
1961: element_props[ek][ej][ei].gp_coords[NSD*p+2] = gp_z;
1962: }
1963: }
1964: }
1965: }
1967: /* define coefficients */
1968: PetscOptionsGetInt(PETSC_NULL,"-model",&model_definition,PETSC_NULL);
1970: switch (model_definition) {
1971: case 0: /* isoviscous */
1972: for (ek = sez; ek < sez+Kmz; ek++) {
1973: for (ej = sey; ej < sey+Jmy; ej++) {
1974: for (ei = sex; ei < sex+Imx; ei++) {
1976: for (p = 0; p < GAUSS_POINTS; p++) {
1977: PetscReal coord_x = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p ]);
1978: PetscReal coord_y = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p+1]);
1979: PetscReal coord_z = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p+2]);
1981: element_props[ek][ej][ei].eta[p] = 1.0;
1983: element_props[ek][ej][ei].fx[p] = 0.0*coord_x;
1984: element_props[ek][ej][ei].fy[p] = -sin((double)2.2*M_PI*coord_y)*cos(1.0*M_PI*coord_x);
1985: element_props[ek][ej][ei].fz[p] = 0.0*coord_z;
1986: element_props[ek][ej][ei].hc[p] = 0.0;
1987: }
1988: }
1989: }
1990: }
1991: break;
1993: case 1: /* manufactured */
1994: for (ek = sez; ek < sez+Kmz; ek++) {
1995: for (ej = sey; ej < sey+Jmy; ej++) {
1996: for (ei = sex; ei < sex+Imx; ei++) {
1997: PetscReal eta,Fm[NSD],Fc,pos[NSD];
1999: for (p = 0; p < GAUSS_POINTS; p++) {
2000: PetscReal coord_x = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p ]);
2001: PetscReal coord_y = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p+1]);
2002: PetscReal coord_z = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p+2]);
2004: pos[0] = coord_x;
2005: pos[1] = coord_y;
2006: pos[2] = coord_z;
2008: evaluate_MS_FrankKamentski(pos,PETSC_NULL,PETSC_NULL,&eta,Fm,&Fc);
2009: element_props[ek][ej][ei].eta[p] = eta;
2011: element_props[ek][ej][ei].fx[p] = Fm[0];
2012: element_props[ek][ej][ei].fy[p] = Fm[1];
2013: element_props[ek][ej][ei].fz[p] = Fm[2];
2014: element_props[ek][ej][ei].hc[p] = Fc;
2015: }
2016: }
2017: }
2018: }
2019: break;
2021: case 2: /* solcx */
2022: for (ek = sez; ek < sez+Kmz; ek++) {
2023: for (ej = sey; ej < sey+Jmy; ej++) {
2024: for (ei = sex; ei < sex+Imx; ei++) {
2025: for (p = 0; p < GAUSS_POINTS; p++) {
2026: PetscReal coord_x = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p ]);
2027: PetscReal coord_y = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p+1]);
2028: PetscReal coord_z = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p+2]);
2030: element_props[ek][ej][ei].eta[p] = 1.0;
2032: element_props[ek][ej][ei].fx[p] = 0.0;
2033: element_props[ek][ej][ei].fy[p] = -sin((double)3*M_PI*coord_y)*cos(1.0*M_PI*coord_x);
2034: element_props[ek][ej][ei].fz[p] = 0.0*coord_z;
2035: element_props[ek][ej][ei].hc[p] = 0.0;
2036: }
2037: }
2038: }
2039: }
2040: break;
2042: case 3: /* sinker */
2043: for (ek = sez; ek < sez+Kmz; ek++) {
2044: for (ej = sey; ej < sey+Jmy; ej++) {
2045: for (ei = sex; ei < sex+Imx; ei++) {
2046: for (p = 0; p < GAUSS_POINTS; p++) {
2047: PetscReal xp = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p ]);
2048: PetscReal yp = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p+1]);
2049: PetscReal zp = PetscRealPart(element_props[ek][ej][ei].gp_coords[NSD*p+2]);
2051: element_props[ek][ej][ei].eta[p] = 1.0e-2;
2052: element_props[ek][ej][ei].fx[p] = 0.0;
2053: element_props[ek][ej][ei].fy[p] = 0.0;
2054: element_props[ek][ej][ei].fz[p] = 0.0;
2055: element_props[ek][ej][ei].hc[p] = 0.0;
2057: if ( (fabs(xp-0.5) < 0.2) &&
2058: (fabs(yp-0.5) < 0.2) &&
2059: (fabs(zp-0.5) < 0.2) ) {
2060: element_props[ek][ej][ei].eta[p] = 1.0;
2061: element_props[ek][ej][ei].fz[p] = 1.0;
2062: }
2064: }
2065: }
2066: }
2067: }
2068: break;
2070: default:
2071: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"No default model is supported. Choose either -model {0,1,2,3}");
2072: break;
2073: }
2075: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
2076: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
2078: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
2079: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
2080: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
2082: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
2083: DMGetMatrix(da_Stokes,MATAIJ,&A);
2084: DMGetMatrix(da_Stokes,MATAIJ,&B);
2085: MatGetVecs(A,&f,&X);
2087: /* assemble A11 */
2088: MatZeroEntries(A);
2089: MatZeroEntries(B);
2090: VecZeroEntries(f);
2092: AssembleA_Stokes(A,da_Stokes,da_prop,properties);
2093: AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
2094: /* build force vector */
2095: AssembleF_Stokes(f,da_Stokes,da_prop,properties);
2097: /* SOLVE */
2098: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
2099: KSPSetOptionsPrefix(ksp_S,"stokes_"); /* stokes */
2100: KSPSetOperators(ksp_S,A,B,SAME_NONZERO_PATTERN);
2101: KSPSetFromOptions(ksp_S);
2103: {
2104: PC pc;
2105: const PetscInt ufields[] = {0,1,2},pfields[] = {3};
2106: KSPGetPC(ksp_S,&pc);
2107: PCFieldSplitSetBlockSize(pc,4);
2108: PCFieldSplitSetFields(pc,"u",3,ufields);
2109: PCFieldSplitSetFields(pc,"p",1,pfields);
2110: }
2112: {
2113: PC pc;
2114: PetscBool same;
2115: KSPGetPC(ksp_S,&pc);
2116: same = PETSC_FALSE;
2117: PetscTypeCompare((PetscObject)pc,PCMG,&same);
2118: if (same) {
2119: PCMGSetupViaCoarsen(pc,da_Stokes);
2120: }
2121: }
2123: {
2124: PetscBool stokes_monitor;
2125: PetscOptionsGetBool(PETSC_NULL,"-stokes_ksp_monitor_blocks",&stokes_monitor,0);
2126: if (stokes_monitor) {
2127: KSPMonitorSet(ksp_S,KSPMonitorStokesBlocks,PETSC_NULL,PETSC_NULL);
2128: }
2129: }
2130: KSPSolve(ksp_S,f,X);
2131: {
2132: PetscBool flg = PETSC_TRUE;
2133: PetscOptionsGetBool(PETSC_NULL,"-write_pvts",&flg,PETSC_NULL);
2134: if (flg) {DAView3DPVTS(da_Stokes,X,"up");}
2135: }
2136: {
2137: PetscBool flg = PETSC_FALSE;
2138: char filename[PETSC_MAX_PATH_LEN];
2139: PetscOptionsGetString(PETSC_NULL,"-write_binary",filename,sizeof filename,&flg);
2140: if (flg) {
2141: PetscViewer viewer;
2142: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer);
2143: VecView(X,viewer);
2144: PetscViewerDestroy(&viewer);
2145: }
2146: }
2147: KSPGetIterationNumber(ksp_S,&its);
2149: /* verify */
2150: if (model_definition == 1) {
2151: DM da_Stokes_analytic;
2152: Vec X_analytic;
2154: DMDACreateManufacturedSolution(mx,my,mz,&da_Stokes_analytic,&X_analytic);
2155: DAView3DPVTS(da_Stokes_analytic,X_analytic,"ms");
2156: DMDAIntegrateErrors3D(da_Stokes_analytic,X,X_analytic);
2157: DAView3DPVTS(da_Stokes,X,"up2");
2159: DMDestroy(&da_Stokes_analytic);
2160: VecDestroy(&X_analytic);
2161: }
2164: KSPDestroy(&ksp_S);
2165: VecDestroy(&X);
2166: VecDestroy(&f);
2167: MatDestroy(&A);
2168: MatDestroy(&B);
2170: DMDestroy(&da_Stokes);
2171: DMDestroy(&da_prop);
2173: VecDestroy(&properties);
2174: VecDestroy(&l_properties);
2176: return(0);
2177: }
2181: int main(int argc,char **args)
2182: {
2184: PetscInt mx,my,mz;
2186: PetscInitialize(&argc,&args,(char *)0,help);
2188: mx = my = mz = 10;
2189: PetscOptionsGetInt(PETSC_NULL,"-mx",&mx,PETSC_NULL);
2190: my = mx; mz = mx;
2191: PetscOptionsGetInt(PETSC_NULL,"-my",&my,PETSC_NULL);
2192: PetscOptionsGetInt(PETSC_NULL,"-mz",&mz,PETSC_NULL);
2194: solve_stokes_3d_coupled(mx,my,mz);
2196: PetscFinalize();
2197: return 0;
2198: }