Actual source code: ex42.c
petsc-3.9.4 2018-09-11
1: static char help[] = "Solves the incompressible, variable viscosity stokes equation in 3d using Q1Q1 elements, \n\
2: stabilized with Bochev's polynomial projection method. Note that implementation here assumes \n\
3: all boundaries are free-slip, i.e. zero normal flow and zero tangential stress \n\
4: -mx : number elements in x-direction \n\
5: -my : number elements in y-direction \n\
6: -mz : number elements in z-direction \n\
7: -stokes_ksp_monitor_blocks : active monitor for each component u,v,w,p \n\
8: -model : defines viscosity and forcing function. Choose either: 0 (isoviscous), 1 (manufactured-broken), 2 (solcx-2d), 3 (sinker), 4 (subdomain jumps) \n";
10: /* Contributed by Dave May */
12: #include <petscksp.h>
13: #include <petscdmda.h>
15: #define PROFILE_TIMING
16: #define ASSEMBLE_LOWER_TRIANGULAR
18: #define NSD 3 /* number of spatial dimensions */
19: #define NODES_PER_EL 8 /* nodes per element */
20: #define U_DOFS 3 /* degrees of freedom per velocity node */
21: #define P_DOFS 1 /* degrees of freedom per pressure node */
22: #define GAUSS_POINTS 8
24: /* Gauss point based evaluation */
25: typedef struct {
26: PetscScalar gp_coords[NSD*GAUSS_POINTS];
27: PetscScalar eta[GAUSS_POINTS];
28: PetscScalar fx[GAUSS_POINTS];
29: PetscScalar fy[GAUSS_POINTS];
30: PetscScalar fz[GAUSS_POINTS];
31: PetscScalar hc[GAUSS_POINTS];
32: } GaussPointCoefficients;
34: typedef struct {
35: PetscScalar u_dof;
36: PetscScalar v_dof;
37: PetscScalar w_dof;
38: PetscScalar p_dof;
39: } StokesDOF;
41: typedef struct _p_CellProperties *CellProperties;
42: struct _p_CellProperties {
43: PetscInt ncells;
44: PetscInt mx,my,mz;
45: PetscInt sex,sey,sez;
46: GaussPointCoefficients *gpc;
47: };
49: /* elements */
50: PetscErrorCode CellPropertiesCreate(DM da_stokes,CellProperties *C)
51: {
53: CellProperties cells;
54: PetscInt mx,my,mz,sex,sey,sez;
57: PetscNew(&cells);
59: DMDAGetElementsCorners(da_stokes,&sex,&sey,&sez);
60: DMDAGetElementsSizes(da_stokes,&mx,&my,&mz);
61: cells->mx = mx;
62: cells->my = my;
63: cells->mz = mz;
64: cells->ncells = mx * my * mz;
65: cells->sex = sex;
66: cells->sey = sey;
67: cells->sez = sez;
69: PetscMalloc1(mx*my*mz,&cells->gpc);
71: *C = cells;
72: return(0);
73: }
75: PetscErrorCode CellPropertiesDestroy(CellProperties *C)
76: {
78: CellProperties cells;
81: if (!C) return(0);
82: cells = *C;
83: PetscFree(cells->gpc);
84: PetscFree(cells);
85: *C = NULL;
86: return(0);
87: }
89: PetscErrorCode CellPropertiesGetCell(CellProperties C,PetscInt II,PetscInt J,PetscInt K,GaussPointCoefficients **G)
90: {
92: *G = &C->gpc[(II-C->sex) + (J-C->sey)*C->mx + (K-C->sez)*C->mx*C->my];
93: return(0);
94: }
96: /* FEM routines */
97: /*
98: Element: Local basis function ordering
99: 1-----2
100: | |
101: | |
102: 0-----3
103: */
104: static void ShapeFunctionQ13D_Evaluate(PetscScalar _xi[],PetscScalar Ni[])
105: {
106: PetscReal xi = PetscRealPart(_xi[0]);
107: PetscReal eta = PetscRealPart(_xi[1]);
108: PetscReal zeta = PetscRealPart(_xi[2]);
110: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
111: Ni[1] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
112: Ni[2] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
113: Ni[3] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
115: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
116: Ni[5] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
117: Ni[6] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
118: Ni[7] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
119: }
121: static void ShapeFunctionQ13D_Evaluate_dxi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
122: {
123: PetscReal xi = PetscRealPart(_xi[0]);
124: PetscReal eta = PetscRealPart(_xi[1]);
125: PetscReal zeta = PetscRealPart(_xi[2]);
126: /* xi deriv */
127: GNi[0][0] = -0.125 * (1.0 - eta) * (1.0 - zeta);
128: GNi[0][1] = -0.125 * (1.0 + eta) * (1.0 - zeta);
129: GNi[0][2] = 0.125 * (1.0 + eta) * (1.0 - zeta);
130: GNi[0][3] = 0.125 * (1.0 - eta) * (1.0 - zeta);
132: GNi[0][4] = -0.125 * (1.0 - eta) * (1.0 + zeta);
133: GNi[0][5] = -0.125 * (1.0 + eta) * (1.0 + zeta);
134: GNi[0][6] = 0.125 * (1.0 + eta) * (1.0 + zeta);
135: GNi[0][7] = 0.125 * (1.0 - eta) * (1.0 + zeta);
136: /* eta deriv */
137: GNi[1][0] = -0.125 * (1.0 - xi) * (1.0 - zeta);
138: GNi[1][1] = 0.125 * (1.0 - xi) * (1.0 - zeta);
139: GNi[1][2] = 0.125 * (1.0 + xi) * (1.0 - zeta);
140: GNi[1][3] = -0.125 * (1.0 + xi) * (1.0 - zeta);
142: GNi[1][4] = -0.125 * (1.0 - xi) * (1.0 + zeta);
143: GNi[1][5] = 0.125 * (1.0 - xi) * (1.0 + zeta);
144: GNi[1][6] = 0.125 * (1.0 + xi) * (1.0 + zeta);
145: GNi[1][7] = -0.125 * (1.0 + xi) * (1.0 + zeta);
146: /* zeta deriv */
147: GNi[2][0] = -0.125 * (1.0 - xi) * (1.0 - eta);
148: GNi[2][1] = -0.125 * (1.0 - xi) * (1.0 + eta);
149: GNi[2][2] = -0.125 * (1.0 + xi) * (1.0 + eta);
150: GNi[2][3] = -0.125 * (1.0 + xi) * (1.0 - eta);
152: GNi[2][4] = 0.125 * (1.0 - xi) * (1.0 - eta);
153: GNi[2][5] = 0.125 * (1.0 - xi) * (1.0 + eta);
154: GNi[2][6] = 0.125 * (1.0 + xi) * (1.0 + eta);
155: GNi[2][7] = 0.125 * (1.0 + xi) * (1.0 - eta);
156: }
158: static void matrix_inverse_3x3(PetscScalar A[3][3],PetscScalar B[3][3])
159: {
160: PetscScalar t4, t6, t8, t10, t12, t14, t17;
162: t4 = A[2][0] * A[0][1];
163: t6 = A[2][0] * A[0][2];
164: t8 = A[1][0] * A[0][1];
165: t10 = A[1][0] * A[0][2];
166: t12 = A[0][0] * A[1][1];
167: t14 = A[0][0] * A[1][2];
168: t17 = 0.1e1 / (t4 * A[1][2] - t6 * A[1][1] - t8 * A[2][2] + t10 * A[2][1] + t12 * A[2][2] - t14 * A[2][1]);
170: B[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) * t17;
171: B[0][1] = -(A[0][1] * A[2][2] - A[0][2] * A[2][1]) * t17;
172: B[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) * t17;
173: B[1][0] = -(-A[2][0] * A[1][2] + A[1][0] * A[2][2]) * t17;
174: B[1][1] = (-t6 + A[0][0] * A[2][2]) * t17;
175: B[1][2] = -(-t10 + t14) * t17;
176: B[2][0] = (-A[2][0] * A[1][1] + A[1][0] * A[2][1]) * t17;
177: B[2][1] = -(-t4 + A[0][0] * A[2][1]) * t17;
178: B[2][2] = (-t8 + t12) * t17;
179: }
181: static void ShapeFunctionQ13D_Evaluate_dx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
182: {
183: PetscScalar J00,J01,J02,J10,J11,J12,J20,J21,J22;
184: PetscInt n;
185: PetscScalar iJ[3][3],JJ[3][3];
187: J00 = J01 = J02 = 0.0;
188: J10 = J11 = J12 = 0.0;
189: J20 = J21 = J22 = 0.0;
190: for (n=0; n<NODES_PER_EL; n++) {
191: PetscScalar cx = coords[NSD*n + 0];
192: PetscScalar cy = coords[NSD*n + 1];
193: PetscScalar cz = coords[NSD*n + 2];
195: /* J_ij = d(x_j) / d(xi_i) */ /* J_ij = \sum _I GNi[j][I} * x_i */
196: J00 = J00 + GNi[0][n] * cx; /* J_xx */
197: J01 = J01 + GNi[0][n] * cy; /* J_xy = dx/deta */
198: J02 = J02 + GNi[0][n] * cz; /* J_xz = dx/dzeta */
200: J10 = J10 + GNi[1][n] * cx; /* J_yx = dy/dxi */
201: J11 = J11 + GNi[1][n] * cy; /* J_yy */
202: J12 = J12 + GNi[1][n] * cz; /* J_yz */
204: J20 = J20 + GNi[2][n] * cx; /* J_zx */
205: J21 = J21 + GNi[2][n] * cy; /* J_zy */
206: J22 = J22 + GNi[2][n] * cz; /* J_zz */
207: }
209: JJ[0][0] = J00; JJ[0][1] = J01; JJ[0][2] = J02;
210: JJ[1][0] = J10; JJ[1][1] = J11; JJ[1][2] = J12;
211: JJ[2][0] = J20; JJ[2][1] = J21; JJ[2][2] = J22;
213: matrix_inverse_3x3(JJ,iJ);
215: *det_J = J00*J11*J22 - J00*J12*J21 - J10*J01*J22 + J10*J02*J21 + J20*J01*J12 - J20*J02*J11;
217: for (n=0; n<NODES_PER_EL; n++) {
218: GNx[0][n] = GNi[0][n]*iJ[0][0] + GNi[1][n]*iJ[0][1] + GNi[2][n]*iJ[0][2];
219: GNx[1][n] = GNi[0][n]*iJ[1][0] + GNi[1][n]*iJ[1][1] + GNi[2][n]*iJ[1][2];
220: GNx[2][n] = GNi[0][n]*iJ[2][0] + GNi[1][n]*iJ[2][1] + GNi[2][n]*iJ[2][2];
221: }
222: }
224: static void ConstructGaussQuadrature3D(PetscInt *ngp,PetscScalar gp_xi[][NSD],PetscScalar gp_weight[])
225: {
226: *ngp = 8;
227: gp_xi[0][0] = -0.57735026919; gp_xi[0][1] = -0.57735026919; gp_xi[0][2] = -0.57735026919;
228: gp_xi[1][0] = -0.57735026919; gp_xi[1][1] = 0.57735026919; gp_xi[1][2] = -0.57735026919;
229: gp_xi[2][0] = 0.57735026919; gp_xi[2][1] = 0.57735026919; gp_xi[2][2] = -0.57735026919;
230: gp_xi[3][0] = 0.57735026919; gp_xi[3][1] = -0.57735026919; gp_xi[3][2] = -0.57735026919;
232: gp_xi[4][0] = -0.57735026919; gp_xi[4][1] = -0.57735026919; gp_xi[4][2] = 0.57735026919;
233: gp_xi[5][0] = -0.57735026919; gp_xi[5][1] = 0.57735026919; gp_xi[5][2] = 0.57735026919;
234: gp_xi[6][0] = 0.57735026919; gp_xi[6][1] = 0.57735026919; gp_xi[6][2] = 0.57735026919;
235: gp_xi[7][0] = 0.57735026919; gp_xi[7][1] = -0.57735026919; gp_xi[7][2] = 0.57735026919;
237: gp_weight[0] = 1.0;
238: gp_weight[1] = 1.0;
239: gp_weight[2] = 1.0;
240: gp_weight[3] = 1.0;
242: gp_weight[4] = 1.0;
243: gp_weight[5] = 1.0;
244: gp_weight[6] = 1.0;
245: gp_weight[7] = 1.0;
246: }
248: /*
249: i,j are the element indices
250: The unknown is a vector quantity.
251: The s[].c is used to indicate the degree of freedom.
252: */
253: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j,PetscInt k)
254: {
255: PetscInt n;
258: /* velocity */
259: n = 0;
260: /* node 0 */
261: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++; /* Vx0 */
262: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++; /* Vy0 */
263: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++; /* Vz0 */
265: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
266: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
267: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
269: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
270: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
271: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
273: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++;
274: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++;
275: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++;
277: /* */
278: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++; /* Vx4 */
279: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++; /* Vy4 */
280: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++; /* Vz4 */
282: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
283: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
284: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
286: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
287: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
288: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
290: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++;
291: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++;
292: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++;
294: /* pressure */
295: n = 0;
297: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++; /* P0 */
298: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
299: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
300: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++;
302: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++; /* P0 */
303: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
304: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
305: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++;
306: return(0);
307: }
309: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords,PetscInt i,PetscInt j,PetscInt k,PetscScalar el_coord[])
310: {
312: /* get coords for the element */
313: el_coord[0] = coords[k][j][i].x;
314: el_coord[1] = coords[k][j][i].y;
315: el_coord[2] = coords[k][j][i].z;
317: el_coord[3] = coords[k][j+1][i].x;
318: el_coord[4] = coords[k][j+1][i].y;
319: el_coord[5] = coords[k][j+1][i].z;
321: el_coord[6] = coords[k][j+1][i+1].x;
322: el_coord[7] = coords[k][j+1][i+1].y;
323: el_coord[8] = coords[k][j+1][i+1].z;
325: el_coord[9] = coords[k][j][i+1].x;
326: el_coord[10] = coords[k][j][i+1].y;
327: el_coord[11] = coords[k][j][i+1].z;
329: el_coord[12] = coords[k+1][j][i].x;
330: el_coord[13] = coords[k+1][j][i].y;
331: el_coord[14] = coords[k+1][j][i].z;
333: el_coord[15] = coords[k+1][j+1][i].x;
334: el_coord[16] = coords[k+1][j+1][i].y;
335: el_coord[17] = coords[k+1][j+1][i].z;
337: el_coord[18] = coords[k+1][j+1][i+1].x;
338: el_coord[19] = coords[k+1][j+1][i+1].y;
339: el_coord[20] = coords[k+1][j+1][i+1].z;
341: el_coord[21] = coords[k+1][j][i+1].x;
342: el_coord[22] = coords[k+1][j][i+1].y;
343: el_coord[23] = coords[k+1][j][i+1].z;
344: return(0);
345: }
347: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field,PetscInt i,PetscInt j,PetscInt k,StokesDOF nodal_fields[])
348: {
350: /* get the nodal fields for u */
351: nodal_fields[0].u_dof = field[k][j][i].u_dof;
352: nodal_fields[0].v_dof = field[k][j][i].v_dof;
353: nodal_fields[0].w_dof = field[k][j][i].w_dof;
355: nodal_fields[1].u_dof = field[k][j+1][i].u_dof;
356: nodal_fields[1].v_dof = field[k][j+1][i].v_dof;
357: nodal_fields[1].w_dof = field[k][j+1][i].w_dof;
359: nodal_fields[2].u_dof = field[k][j+1][i+1].u_dof;
360: nodal_fields[2].v_dof = field[k][j+1][i+1].v_dof;
361: nodal_fields[2].w_dof = field[k][j+1][i+1].w_dof;
363: nodal_fields[3].u_dof = field[k][j][i+1].u_dof;
364: nodal_fields[3].v_dof = field[k][j][i+1].v_dof;
365: nodal_fields[3].w_dof = field[k][j][i+1].w_dof;
367: nodal_fields[4].u_dof = field[k+1][j][i].u_dof;
368: nodal_fields[4].v_dof = field[k+1][j][i].v_dof;
369: nodal_fields[4].w_dof = field[k+1][j][i].w_dof;
371: nodal_fields[5].u_dof = field[k+1][j+1][i].u_dof;
372: nodal_fields[5].v_dof = field[k+1][j+1][i].v_dof;
373: nodal_fields[5].w_dof = field[k+1][j+1][i].w_dof;
375: nodal_fields[6].u_dof = field[k+1][j+1][i+1].u_dof;
376: nodal_fields[6].v_dof = field[k+1][j+1][i+1].v_dof;
377: nodal_fields[6].w_dof = field[k+1][j+1][i+1].w_dof;
379: nodal_fields[7].u_dof = field[k+1][j][i+1].u_dof;
380: nodal_fields[7].v_dof = field[k+1][j][i+1].v_dof;
381: nodal_fields[7].w_dof = field[k+1][j][i+1].w_dof;
383: /* get the nodal fields for p */
384: nodal_fields[0].p_dof = field[k][j][i].p_dof;
385: nodal_fields[1].p_dof = field[k][j+1][i].p_dof;
386: nodal_fields[2].p_dof = field[k][j+1][i+1].p_dof;
387: nodal_fields[3].p_dof = field[k][j][i+1].p_dof;
389: nodal_fields[4].p_dof = field[k+1][j][i].p_dof;
390: nodal_fields[5].p_dof = field[k+1][j+1][i].p_dof;
391: nodal_fields[6].p_dof = field[k+1][j+1][i+1].p_dof;
392: nodal_fields[7].p_dof = field[k+1][j][i+1].p_dof;
393: return(0);
394: }
396: static PetscInt ASS_MAP_wIwDI_uJuDJ(PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)
397: {
398: PetscInt ij;
399: PETSC_UNUSED PetscInt r,c,nr,nc;
401: nr = w_NPE*w_dof;
402: nc = u_NPE*u_dof;
404: r = w_dof*wi+wd;
405: c = u_dof*ui+ud;
407: ij = r*nc+c;
409: return ij;
410: }
412: static PetscErrorCode DMDASetValuesLocalStencil3D_ADD_VALUES(StokesDOF ***fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
413: {
414: PetscInt n,II,J,K;
417: for (n = 0; n<NODES_PER_EL; n++) {
418: II = u_eqn[NSD*n].i;
419: J = u_eqn[NSD*n].j;
420: K = u_eqn[NSD*n].k;
422: fields_F[K][J][II].u_dof = fields_F[K][J][II].u_dof+Fe_u[NSD*n];
424: II = u_eqn[NSD*n+1].i;
425: J = u_eqn[NSD*n+1].j;
426: K = u_eqn[NSD*n+1].k;
428: fields_F[K][J][II].v_dof = fields_F[K][J][II].v_dof+Fe_u[NSD*n+1];
430: II = u_eqn[NSD*n+2].i;
431: J = u_eqn[NSD*n+2].j;
432: K = u_eqn[NSD*n+2].k;
433: fields_F[K][J][II].w_dof = fields_F[K][J][II].w_dof+Fe_u[NSD*n+2];
435: II = p_eqn[n].i;
436: J = p_eqn[n].j;
437: K = p_eqn[n].k;
439: fields_F[K][J][II].p_dof = fields_F[K][J][II].p_dof+Fe_p[n];
441: }
442: return(0);
443: }
445: static void FormStressOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
446: {
447: PetscInt ngp;
448: PetscScalar gp_xi[GAUSS_POINTS][NSD];
449: PetscScalar gp_weight[GAUSS_POINTS];
450: PetscInt p,i,j,k;
451: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
452: PetscScalar J_p,tildeD[6];
453: PetscScalar B[6][U_DOFS*NODES_PER_EL];
454: const PetscInt nvdof = U_DOFS*NODES_PER_EL;
456: /* define quadrature rule */
457: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
459: /* evaluate integral */
460: for (p = 0; p < ngp; p++) {
461: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
462: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
464: for (i = 0; i < NODES_PER_EL; i++) {
465: PetscScalar d_dx_i = GNx_p[0][i];
466: PetscScalar d_dy_i = GNx_p[1][i];
467: PetscScalar d_dz_i = GNx_p[2][i];
469: B[0][3*i] = d_dx_i; B[0][3*i+1] = 0.0; B[0][3*i+2] = 0.0;
470: B[1][3*i] = 0.0; B[1][3*i+1] = d_dy_i; B[1][3*i+2] = 0.0;
471: B[2][3*i] = 0.0; B[2][3*i+1] = 0.0; B[2][3*i+2] = d_dz_i;
473: B[3][3*i] = d_dy_i; B[3][3*i+1] = d_dx_i; B[3][3*i+2] = 0.0; /* e_xy */
474: B[4][3*i] = d_dz_i; B[4][3*i+1] = 0.0; B[4][3*i+2] = d_dx_i; /* e_xz */
475: B[5][3*i] = 0.0; B[5][3*i+1] = d_dz_i; B[5][3*i+2] = d_dy_i; /* e_yz */
476: }
479: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
480: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
481: tildeD[2] = 2.0*gp_weight[p]*J_p*eta[p];
483: tildeD[3] = gp_weight[p]*J_p*eta[p];
484: tildeD[4] = gp_weight[p]*J_p*eta[p];
485: tildeD[5] = gp_weight[p]*J_p*eta[p];
487: /* form Bt tildeD B */
488: /*
489: Ke_ij = Bt_ik . D_kl . B_lj
490: = B_ki . D_kl . B_lj
491: = B_ki . D_kk . B_kj
492: */
493: for (i = 0; i < nvdof; i++) {
494: for (j = i; j < nvdof; j++) {
495: for (k = 0; k < 6; k++) {
496: Ke[i*nvdof+j] += B[k][i]*tildeD[k]*B[k][j];
497: }
498: }
499: }
501: }
502: /* fill lower triangular part */
503: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
504: for (i = 0; i < nvdof; i++) {
505: for (j = i; j < nvdof; j++) {
506: Ke[j*nvdof+i] = Ke[i*nvdof+j];
507: }
508: }
509: #endif
510: }
512: static void FormGradientOperatorQ13D(PetscScalar Ke[],PetscScalar coords[])
513: {
514: PetscInt ngp;
515: PetscScalar gp_xi[GAUSS_POINTS][NSD];
516: PetscScalar gp_weight[GAUSS_POINTS];
517: PetscInt p,i,j,di;
518: PetscScalar Ni_p[NODES_PER_EL];
519: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
520: PetscScalar J_p,fac;
522: /* define quadrature rule */
523: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
525: /* evaluate integral */
526: for (p = 0; p < ngp; p++) {
527: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
528: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
529: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
530: fac = gp_weight[p]*J_p;
532: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
533: for (di = 0; di < NSD; di++) { /* u dofs */
534: for (j = 0; j < NODES_PER_EL; j++) { /* p nodes, p dofs = 1 (ie no loop) */
535: PetscInt IJ;
536: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,3,j,0,NODES_PER_EL,1);
538: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
539: }
540: }
541: }
542: }
543: }
545: static void FormDivergenceOperatorQ13D(PetscScalar De[],PetscScalar coords[])
546: {
547: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
548: PetscInt i,j;
549: PetscInt nr_g,nc_g;
551: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
552: FormGradientOperatorQ13D(Ge,coords);
554: nr_g = U_DOFS*NODES_PER_EL;
555: nc_g = P_DOFS*NODES_PER_EL;
557: for (i = 0; i < nr_g; i++) {
558: for (j = 0; j < nc_g; j++) {
559: De[nr_g*j+i] = Ge[nc_g*i+j];
560: }
561: }
562: }
564: static void FormStabilisationOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
565: {
566: PetscInt ngp;
567: PetscScalar gp_xi[GAUSS_POINTS][NSD];
568: PetscScalar gp_weight[GAUSS_POINTS];
569: PetscInt p,i,j;
570: PetscScalar Ni_p[NODES_PER_EL];
571: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
572: PetscScalar J_p,fac,eta_avg;
574: /* define quadrature rule */
575: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
577: /* evaluate integral */
578: for (p = 0; p < ngp; p++) {
579: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
580: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
581: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
582: fac = gp_weight[p]*J_p;
583: /*
584: for (i = 0; i < NODES_PER_EL; i++) {
585: for (j = i; j < NODES_PER_EL; j++) {
586: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
587: }
588: }
589: */
591: for (i = 0; i < NODES_PER_EL; i++) {
592: for (j = 0; j < NODES_PER_EL; j++) {
593: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
594: }
595: }
596: }
598: /* scale */
599: eta_avg = 0.0;
600: for (p = 0; p < ngp; p++) eta_avg += eta[p];
601: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
602: fac = 1.0/eta_avg;
604: /*
605: for (i = 0; i < NODES_PER_EL; i++) {
606: for (j = i; j < NODES_PER_EL; j++) {
607: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
608: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
609: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
610: #endif
611: }
612: }
613: */
615: for (i = 0; i < NODES_PER_EL; i++) {
616: for (j = 0; j < NODES_PER_EL; j++) {
617: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
618: }
619: }
620: }
622: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
623: {
624: PetscInt ngp;
625: PetscScalar gp_xi[GAUSS_POINTS][NSD];
626: PetscScalar gp_weight[GAUSS_POINTS];
627: PetscInt p,i,j;
628: PetscScalar Ni_p[NODES_PER_EL];
629: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
630: PetscScalar J_p,fac,eta_avg;
632: /* define quadrature rule */
633: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
635: /* evaluate integral */
636: for (p = 0; p < ngp; p++) {
637: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
638: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
639: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
640: fac = gp_weight[p]*J_p;
642: /*
643: for (i = 0; i < NODES_PER_EL; i++) {
644: for (j = i; j < NODES_PER_EL; j++) {
645: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
646: }
647: }
648: */
650: for (i = 0; i < NODES_PER_EL; i++) {
651: for (j = 0; j < NODES_PER_EL; j++) {
652: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
653: }
654: }
655: }
657: /* scale */
658: eta_avg = 0.0;
659: for (p = 0; p < ngp; p++) eta_avg += eta[p];
660: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
661: fac = 1.0/eta_avg;
662: /*
663: for (i = 0; i < NODES_PER_EL; i++) {
664: for (j = i; j < NODES_PER_EL; j++) {
665: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
666: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
667: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
668: #endif
669: }
670: }
671: */
673: for (i = 0; i < NODES_PER_EL; i++) {
674: for (j = 0; j < NODES_PER_EL; j++) {
675: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
676: }
677: }
678: }
680: static void FormMomentumRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[],PetscScalar fz[])
681: {
682: PetscInt ngp;
683: PetscScalar gp_xi[GAUSS_POINTS][NSD];
684: PetscScalar gp_weight[GAUSS_POINTS];
685: PetscInt p,i;
686: PetscScalar Ni_p[NODES_PER_EL];
687: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
688: PetscScalar J_p,fac;
690: /* define quadrature rule */
691: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
693: /* evaluate integral */
694: for (p = 0; p < ngp; p++) {
695: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
696: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
697: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
698: fac = gp_weight[p]*J_p;
700: for (i = 0; i < NODES_PER_EL; i++) {
701: Fe[NSD*i] -= fac*Ni_p[i]*fx[p];
702: Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
703: Fe[NSD*i+2] -= fac*Ni_p[i]*fz[p];
704: }
705: }
706: }
708: static void FormContinuityRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar hc[])
709: {
710: PetscInt ngp;
711: PetscScalar gp_xi[GAUSS_POINTS][NSD];
712: PetscScalar gp_weight[GAUSS_POINTS];
713: PetscInt p,i;
714: PetscScalar Ni_p[NODES_PER_EL];
715: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
716: PetscScalar J_p,fac;
718: /* define quadrature rule */
719: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
721: /* evaluate integral */
722: for (p = 0; p < ngp; p++) {
723: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
724: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
725: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
726: fac = gp_weight[p]*J_p;
728: for (i = 0; i < NODES_PER_EL; i++) Fe[i] -= fac*Ni_p[i]*hc[p];
729: }
730: }
732: #define _ZERO_ROWCOL_i(A,i) { \
733: PetscInt KK; \
734: PetscScalar tmp = A[24*(i)+(i)]; \
735: for (KK=0;KK<24;KK++) A[24*(i)+KK]=0.0; \
736: for (KK=0;KK<24;KK++) A[24*KK+(i)]=0.0; \
737: A[24*(i)+(i)] = tmp;} \
739: #define _ZERO_ROW_i(A,i) { \
740: PetscInt KK; \
741: for (KK=0;KK<8;KK++) A[8*(i)+KK]=0.0;}
743: #define _ZERO_COL_i(A,i) { \
744: PetscInt KK; \
745: for (KK=0;KK<8;KK++) A[24*KK+(i)]=0.0;}
747: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,CellProperties cell_properties)
748: {
749: DM cda;
750: Vec coords;
751: DMDACoor3d ***_coords;
752: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
753: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
754: PetscInt sex,sey,sez,mx,my,mz;
755: PetscInt ei,ej,ek;
756: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
757: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
758: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
759: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
760: PetscScalar el_coords[NODES_PER_EL*NSD];
761: GaussPointCoefficients *props;
762: PetscScalar *prop_eta;
763: PetscInt n,M,N,P;
764: PetscErrorCode ierr;
767: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
768: /* setup for coords */
769: DMGetCoordinateDM(stokes_da,&cda);
770: DMGetCoordinatesLocal(stokes_da,&coords);
771: DMDAVecGetArray(cda,coords,&_coords);
772: DMDAGetElementsCorners(stokes_da,&sex,&sey,&sez);
773: DMDAGetElementsSizes(stokes_da,&mx,&my,&mz);
774: for (ek = sez; ek < sez+mz; ek++) {
775: for (ej = sey; ej < sey+my; ej++) {
776: for (ei = sex; ei < sex+mx; ei++) {
777: /* get coords for the element */
778: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
779: /* get cell properties */
780: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
781: /* get coefficients for the element */
782: prop_eta = props->eta;
784: /* initialise element stiffness matrix */
785: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
786: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
787: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
788: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
790: /* form element stiffness matrix */
791: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
792: FormGradientOperatorQ13D(Ge,el_coords);
793: /*#if defined(ASSEMBLE_LOWER_TRIANGULAR)*/
794: FormDivergenceOperatorQ13D(De,el_coords);
795: /*#endif*/
796: FormStabilisationOperatorQ13D(Ce,el_coords,prop_eta);
798: /* insert element matrix into global matrix */
799: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
801: for (n=0; n<NODES_PER_EL; n++) {
802: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
803: _ZERO_ROWCOL_i(Ae,3*n);
804: _ZERO_ROW_i (Ge,3*n);
805: _ZERO_COL_i (De,3*n);
806: }
808: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
809: _ZERO_ROWCOL_i(Ae,3*n+1);
810: _ZERO_ROW_i (Ge,3*n+1);
811: _ZERO_COL_i (De,3*n+1);
812: }
814: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
815: _ZERO_ROWCOL_i(Ae,3*n+2);
816: _ZERO_ROW_i (Ge,3*n+2);
817: _ZERO_COL_i (De,3*n+2);
818: }
819: }
820: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
821: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
822: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
823: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
824: }
825: }
826: }
827: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
828: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
830: DMDAVecRestoreArray(cda,coords,&_coords);
832: return(0);
833: }
835: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,CellProperties cell_properties)
836: {
837: DM cda;
838: Vec coords;
839: DMDACoor3d ***_coords;
840: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
841: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
842: PetscInt sex,sey,sez,mx,my,mz;
843: PetscInt ei,ej,ek;
844: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
845: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
846: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
847: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
848: PetscScalar el_coords[NODES_PER_EL*NSD];
849: GaussPointCoefficients *props;
850: PetscScalar *prop_eta;
851: PetscInt n,M,N,P;
852: PetscErrorCode ierr;
855: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
856: /* setup for coords */
857: DMGetCoordinateDM(stokes_da,&cda);
858: DMGetCoordinatesLocal(stokes_da,&coords);
859: DMDAVecGetArray(cda,coords,&_coords);
861: DMDAGetElementsCorners(stokes_da,&sex,&sey,&sez);
862: DMDAGetElementsSizes(stokes_da,&mx,&my,&mz);
863: for (ek = sez; ek < sez+mz; ek++) {
864: for (ej = sey; ej < sey+my; ej++) {
865: for (ei = sex; ei < sex+mx; ei++) {
866: /* get coords for the element */
867: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
868: /* get cell properties */
869: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
870: /* get coefficients for the element */
871: prop_eta = props->eta;
873: /* initialise element stiffness matrix */
874: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
875: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
876: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
877: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
879: /* form element stiffness matrix */
880: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
881: FormGradientOperatorQ13D(Ge,el_coords);
882: /* FormDivergenceOperatorQ13D(De,el_coords); */
883: FormScaledMassMatrixOperatorQ13D(Ce,el_coords,prop_eta);
885: /* insert element matrix into global matrix */
886: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
888: for (n=0; n<NODES_PER_EL; n++) {
889: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
890: _ZERO_ROWCOL_i(Ae,3*n);
891: _ZERO_ROW_i (Ge,3*n);
892: _ZERO_COL_i (De,3*n);
893: }
895: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
896: _ZERO_ROWCOL_i(Ae,3*n+1);
897: _ZERO_ROW_i (Ge,3*n+1);
898: _ZERO_COL_i (De,3*n+1);
899: }
901: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
902: _ZERO_ROWCOL_i(Ae,3*n+2);
903: _ZERO_ROW_i (Ge,3*n+2);
904: _ZERO_COL_i (De,3*n+2);
905: }
906: }
907: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
908: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
909: /*MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);*/
910: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
911: }
912: }
913: }
914: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
915: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
917: DMDAVecRestoreArray(cda,coords,&_coords);
918: return(0);
919: }
921: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,CellProperties cell_properties)
922: {
923: DM cda;
924: Vec coords;
925: DMDACoor3d ***_coords;
926: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
927: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
928: PetscInt sex,sey,sez,mx,my,mz;
929: PetscInt ei,ej,ek;
930: PetscScalar Fe[NODES_PER_EL*U_DOFS];
931: PetscScalar He[NODES_PER_EL*P_DOFS];
932: PetscScalar el_coords[NODES_PER_EL*NSD];
933: GaussPointCoefficients *props;
934: PetscScalar *prop_fx,*prop_fy,*prop_fz,*prop_hc;
935: Vec local_F;
936: StokesDOF ***ff;
937: PetscInt n,M,N,P;
938: PetscErrorCode ierr;
941: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
942: /* setup for coords */
943: DMGetCoordinateDM(stokes_da,&cda);
944: DMGetCoordinatesLocal(stokes_da,&coords);
945: DMDAVecGetArray(cda,coords,&_coords);
947: /* get acces to the vector */
948: DMGetLocalVector(stokes_da,&local_F);
949: VecZeroEntries(local_F);
950: DMDAVecGetArray(stokes_da,local_F,&ff);
951: DMDAGetElementsCorners(stokes_da,&sex,&sey,&sez);
952: DMDAGetElementsSizes(stokes_da,&mx,&my,&mz);
953: for (ek = sez; ek < sez+mz; ek++) {
954: for (ej = sey; ej < sey+my; ej++) {
955: for (ei = sex; ei < sex+mx; ei++) {
956: /* get coords for the element */
957: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
958: /* get cell properties */
959: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
960: /* get coefficients for the element */
961: prop_fx = props->fx;
962: prop_fy = props->fy;
963: prop_fz = props->fz;
964: prop_hc = props->hc;
966: /* initialise element stiffness matrix */
967: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
968: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
970: /* form element stiffness matrix */
971: FormMomentumRhsQ13D(Fe,el_coords,prop_fx,prop_fy,prop_fz);
972: FormContinuityRhsQ13D(He,el_coords,prop_hc);
974: /* insert element matrix into global matrix */
975: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
977: for (n=0; n<NODES_PER_EL; n++) {
978: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) Fe[3*n] = 0.0;
980: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) Fe[3*n+1] = 0.0;
982: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) Fe[3*n+2] = 0.0;
983: }
985: DMDASetValuesLocalStencil3D_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
986: }
987: }
988: }
989: DMDAVecRestoreArray(stokes_da,local_F,&ff);
990: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
991: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
992: DMRestoreLocalVector(stokes_da,&local_F);
994: DMDAVecRestoreArray(cda,coords,&_coords);
995: return(0);
996: }
998: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta,PetscReal *MX,PetscReal *MY,PetscReal *MZ)
999: {
1000: *theta = 0.0;
1001: *MX = 2.0 * PETSC_PI;
1002: *MY = 2.0 * PETSC_PI;
1003: *MZ = 2.0 * PETSC_PI;
1004: }
1005: static void evaluate_MS_FrankKamentski(PetscReal pos[],PetscReal v[],PetscReal *p,PetscReal *eta,PetscReal Fm[],PetscReal *Fc)
1006: {
1007: PetscReal x,y,z;
1008: PetscReal theta,MX,MY,MZ;
1010: evaluate_MS_FrankKamentski_constants(&theta,&MX,&MY,&MZ);
1011: x = pos[0];
1012: y = pos[1];
1013: z = pos[2];
1014: if (v) {
1015: /*
1016: v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1017: v[1] = z*cos(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1018: v[2] = -(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z);
1019: */
1020: /*
1021: v[0] = PetscSinReal(PETSC_PI*x);
1022: v[1] = PetscSinReal(PETSC_PI*y);
1023: v[2] = PetscSinReal(PETSC_PI*z);
1024: */
1025: v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1026: v[1] = z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1027: v[2] = PetscPowRealInt(z,2)*(PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)/2 - PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)/2) - PETSC_PI*PetscPowRealInt(z,4)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y)/4;
1028: }
1029: if (p) *p = PetscPowRealInt(x,2) + PetscPowRealInt(y,2) + PetscPowRealInt(z,2);
1030: if (eta) {
1031: /**eta = PetscExpReal(-theta*(1.0 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));*/
1032: *eta = 1.0;
1033: }
1034: if (Fm) {
1035: /*
1036: Fm[0] = -2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(PETSC_PI*x) - 0.2*MZ*theta*(-1.5*PetscPowRealInt(x,2)*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PetscPowRealInt(z,2)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x))*PetscCosReal(MX*x)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) - 0.2*PETSC_PI*MX*theta*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscCosReal(MZ*z)*PetscExpReal(y)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) + 2.0*(3.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 3.0*PETSC_PI*PetscPowRealInt(x,2)*PetscCosReal(2.0*PETSC_PI*z))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 1.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) ;
1037: Fm[1] = -2*y - 0.2*MX*theta*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x))*PetscCosReal(MZ*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) - 0.2*MZ*theta*(-1.5*PetscPowRealInt(y,2)*PetscSinReal(2.0*PETSC_PI*z) + 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y))*PetscCosReal(MX*x)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) + 2.0*(-2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.5*PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 2*PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(-z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) - 6.0*PETSC_PI*PetscPowRealInt(y,2)*PetscCosReal(2.0*PETSC_PI*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));
1038: Fm[2] = -2*z + 8.0*PetscPowRealInt(PETSC_PI,2)*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(2.0*PETSC_PI*z) - 0.2*MX*theta*(-1.5*PetscPowRealInt(x,2)*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PetscPowRealInt(z,2)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x))*PetscCosReal(MZ*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) + 0.4*PETSC_PI*MZ*theta*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(MX*x)*PetscCosReal(2.0*PETSC_PI*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) + 2.0*(-3.0*x*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(-3.0*y*PetscSinReal(2.0*PETSC_PI*z) - 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.5*PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(-1.5*PetscPowRealInt(y,2)*PetscSinReal(2.0*PETSC_PI*z) + 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) ;
1039: */
1040: /*
1041: Fm[0]=-2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*x);
1042: Fm[1]=-2*y - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*y);
1043: Fm[2]=-2*z - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*z);
1044: */
1045: /*
1046: Fm[0] = -2*x + PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 6.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 6.0*PETSC_PI*PetscPowRealInt(x,2)*PetscCosReal(2.0*PETSC_PI*z) - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 2.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x) ;
1047: Fm[1] = -2*y - 6.0*PETSC_PI*PetscPowRealInt(y,2)*PetscCosReal(2.0*PETSC_PI*z) + 2.0*z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 6.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) - 4.0*PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);
1048: Fm[2] = -2*z - 6.0*x*PetscSinReal(2.0*PETSC_PI*z) - 6.0*y*PetscSinReal(2.0*PETSC_PI*z) - PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 8.0*PetscPowRealInt(PETSC_PI,2)*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z) + 3.0*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y) ;
1049: */
1050: Fm[0] = -2*x + 2*z*(PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x) - 1.0*PETSC_PI*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x)) + PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 6.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 2.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x);
1051: Fm[1] = -2*y + 2*z*(-PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)/2 + PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)/2 + PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)) + 2.0*z*PetscCosReal(2.0*PETSC_PI *x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 6.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 4.0*PETSC_PI*z*PetscCosReal(PETSC_PI *y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);
1052: Fm[2] = -2*z + PetscPowRealInt(z,2)*(-2.0*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 2.0*PetscPowRealInt(PETSC_PI,3)*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)) + PetscPowRealInt(z,2)*(PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)/2 - 3*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)/2 + PetscPowRealInt(PETSC_PI,3)*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)/2 - 3*PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)/2) + 1.0*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.25*PetscPowRealInt(PETSC_PI,3)*PetscPowRealInt(z,4)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) - 0.25*PETSC_PI*PetscPowRealInt(z,4)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) - 3.0*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) - 1.0*PETSC_PI*PetscCosReal(PETSC_PI *y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);
1053: }
1054: if (Fc) {
1055: /**Fc = -2.0*PETSC_PI*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(2.0*PETSC_PI*z) - z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y) ;*/
1056: /**Fc = PETSC_PI*PetscCosReal(PETSC_PI*x) + PETSC_PI*PetscCosReal(PETSC_PI*y) + PETSC_PI*PetscCosReal(PETSC_PI*z);*/
1057: /**Fc = -2.0*PETSC_PI*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(2.0*PETSC_PI*z) - z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);*/
1058: *Fc = 0.0;
1059: }
1060: }
1062: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx,PetscInt my,PetscInt mz,DM *_da,Vec *_X)
1063: {
1064: DM da,cda;
1065: Vec X;
1066: StokesDOF ***_stokes;
1067: Vec coords;
1068: DMDACoor3d ***_coords;
1069: PetscInt si,sj,sk,ei,ej,ek,i,j,k;
1073: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1074: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,4,1,NULL,NULL,NULL,&da);
1075: DMSetFromOptions(da);
1076: DMSetUp(da);
1077: DMDASetFieldName(da,0,"anlytic_Vx");
1078: DMDASetFieldName(da,1,"anlytic_Vy");
1079: DMDASetFieldName(da,2,"anlytic_Vz");
1080: DMDASetFieldName(da,3,"analytic_P");
1082: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
1084: DMGetCoordinatesLocal(da,&coords);
1085: DMGetCoordinateDM(da,&cda);
1086: DMDAVecGetArray(cda,coords,&_coords);
1088: DMCreateGlobalVector(da,&X);
1089: DMDAVecGetArray(da,X,&_stokes);
1091: DMDAGetCorners(da,&si,&sj,&sk,&ei,&ej,&ek);
1092: for (k = sk; k < sk+ek; k++) {
1093: for (j = sj; j < sj+ej; j++) {
1094: for (i = si; i < si+ei; i++) {
1095: PetscReal pos[NSD],pressure,vel[NSD];
1097: pos[0] = PetscRealPart(_coords[k][j][i].x);
1098: pos[1] = PetscRealPart(_coords[k][j][i].y);
1099: pos[2] = PetscRealPart(_coords[k][j][i].z);
1101: evaluate_MS_FrankKamentski(pos,vel,&pressure,NULL,NULL,NULL);
1103: _stokes[k][j][i].u_dof = vel[0];
1104: _stokes[k][j][i].v_dof = vel[1];
1105: _stokes[k][j][i].w_dof = vel[2];
1106: _stokes[k][j][i].p_dof = pressure;
1107: }
1108: }
1109: }
1110: DMDAVecRestoreArray(da,X,&_stokes);
1111: DMDAVecRestoreArray(cda,coords,&_coords);
1113: *_da = da;
1114: *_X = X;
1115: return(0);
1116: }
1118: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da,Vec X,Vec X_analytic)
1119: {
1120: DM cda;
1121: Vec coords,X_analytic_local,X_local;
1122: DMDACoor3d ***_coords;
1123: PetscInt sex,sey,sez,mx,my,mz;
1124: PetscInt ei,ej,ek;
1125: PetscScalar el_coords[NODES_PER_EL*NSD];
1126: StokesDOF ***stokes_analytic,***stokes;
1127: StokesDOF stokes_analytic_e[NODES_PER_EL],stokes_e[NODES_PER_EL];
1129: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1130: PetscScalar Ni_p[NODES_PER_EL];
1131: PetscInt ngp;
1132: PetscScalar gp_xi[GAUSS_POINTS][NSD];
1133: PetscScalar gp_weight[GAUSS_POINTS];
1134: PetscInt p,i;
1135: PetscScalar J_p,fac;
1136: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1137: PetscScalar tint_p_ms,tint_p,int_p_ms,int_p;
1138: PetscInt M;
1139: PetscReal xymin[NSD],xymax[NSD];
1143: /* define quadrature rule */
1144: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1146: /* setup for coords */
1147: DMGetCoordinateDM(stokes_da,&cda);
1148: DMGetCoordinatesLocal(stokes_da,&coords);
1149: DMDAVecGetArray(cda,coords,&_coords);
1151: /* setup for analytic */
1152: DMCreateLocalVector(stokes_da,&X_analytic_local);
1153: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1154: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1155: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1157: /* setup for solution */
1158: DMCreateLocalVector(stokes_da,&X_local);
1159: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1160: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1161: DMDAVecGetArray(stokes_da,X_local,&stokes);
1163: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1164: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1166: h = (xymax[0]-xymin[0])/((PetscReal)(M-1));
1168: tp_L2 = tu_L2 = tu_H1 = 0.0;
1169: tint_p_ms = tint_p = 0.0;
1171: DMDAGetElementsCorners(stokes_da,&sex,&sey,&sez);
1172: DMDAGetElementsSizes(stokes_da,&mx,&my,&mz);
1173: for (ek = sez; ek < sez+mz; ek++) {
1174: for (ej = sey; ej < sey+my; ej++) {
1175: for (ei = sex; ei < sex+mx; ei++) {
1176: /* get coords for the element */
1177: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1178: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1179: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1181: /* evaluate integral */
1182: for (p = 0; p < ngp; p++) {
1183: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1184: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1185: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1186: fac = gp_weight[p]*J_p;
1188: for (i = 0; i < NODES_PER_EL; i++) {
1189: tint_p_ms = tint_p_ms+fac*Ni_p[i]*stokes_analytic_e[i].p_dof;
1190: tint_p = tint_p +fac*Ni_p[i]*stokes_e[i].p_dof;
1191: }
1192: }
1193: }
1194: }
1195: }
1196: MPI_Allreduce(&tint_p_ms,&int_p_ms,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1197: MPI_Allreduce(&tint_p,&int_p,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1199: PetscPrintf(PETSC_COMM_WORLD,"\\int P dv %1.4e (h) %1.4e (ms)\n",PetscRealPart(int_p),PetscRealPart(int_p_ms));
1201: /* remove mine and add manufacture one */
1202: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1203: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1205: {
1206: PetscInt k,L,dof;
1207: PetscScalar *fields;
1208: DMDAGetInfo(stokes_da,0, 0,0,0, 0,0,0, &dof,0,0,0,0,0);
1210: VecGetLocalSize(X_local,&L);
1211: VecGetArray(X_local,&fields);
1212: for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1213: VecRestoreArray(X_local,&fields);
1215: VecGetLocalSize(X,&L);
1216: VecGetArray(X,&fields);
1217: for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1218: VecRestoreArray(X,&fields);
1219: }
1221: DMDAVecGetArray(stokes_da,X_local,&stokes);
1222: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1224: for (ek = sez; ek < sez+mz; ek++) {
1225: for (ej = sey; ej < sey+my; ej++) {
1226: for (ei = sex; ei < sex+mx; ei++) {
1227: /* get coords for the element */
1228: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1229: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1230: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1232: /* evaluate integral */
1233: p_e_L2 = 0.0;
1234: u_e_L2 = 0.0;
1235: u_e_H1 = 0.0;
1236: for (p = 0; p < ngp; p++) {
1237: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1238: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1239: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1240: fac = gp_weight[p]*J_p;
1242: for (i = 0; i < NODES_PER_EL; i++) {
1243: PetscScalar u_error,v_error,w_error;
1245: p_e_L2 = p_e_L2+fac*Ni_p[i]*(stokes_e[i].p_dof-stokes_analytic_e[i].p_dof)*(stokes_e[i].p_dof-stokes_analytic_e[i].p_dof);
1247: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1248: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1249: w_error = stokes_e[i].w_dof-stokes_analytic_e[i].w_dof;
1250: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error+w_error*w_error);
1252: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1253: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error
1254: +GNx_p[2][i]*u_error*GNx_p[2][i]*u_error
1255: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1256: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error
1257: +GNx_p[2][i]*v_error*GNx_p[2][i]*v_error
1258: +GNx_p[0][i]*w_error*GNx_p[0][i]*w_error /* dw/dx */
1259: +GNx_p[1][i]*w_error*GNx_p[1][i]*w_error
1260: +GNx_p[2][i]*w_error*GNx_p[2][i]*w_error);
1261: }
1262: }
1264: tp_L2 += p_e_L2;
1265: tu_L2 += u_e_L2;
1266: tu_H1 += u_e_H1;
1267: }
1268: }
1269: }
1270: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1271: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1272: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1273: p_L2 = PetscSqrtScalar(p_L2);
1274: u_L2 = PetscSqrtScalar(u_L2);
1275: u_H1 = PetscSqrtScalar(u_H1);
1277: PetscPrintf(PETSC_COMM_WORLD,"%1.4e %1.4e %1.4e %1.4e \n",PetscRealPart(h),PetscRealPart(p_L2),PetscRealPart(u_L2),PetscRealPart(u_H1));
1280: DMDAVecRestoreArray(cda,coords,&_coords);
1282: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1283: VecDestroy(&X_analytic_local);
1284: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1285: VecDestroy(&X_local);
1286: return(0);
1287: }
1289: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da,Vec FIELD,const char file_prefix[])
1290: {
1291: char vtk_filename[PETSC_MAX_PATH_LEN];
1292: PetscMPIInt rank;
1293: MPI_Comm comm;
1294: FILE *vtk_fp = NULL;
1295: PetscInt si,sj,sk,nx,ny,nz,i;
1296: PetscInt f,n_fields,N;
1297: DM cda;
1298: Vec coords;
1299: Vec l_FIELD;
1300: PetscScalar *_L_FIELD;
1301: PetscInt memory_offset;
1302: PetscScalar *buffer;
1307: /* create file name */
1308: PetscObjectGetComm((PetscObject)da,&comm);
1309: MPI_Comm_rank(comm,&rank);
1310: PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"subdomain-%s-p%1.4d.vts",file_prefix,rank);
1312: /* open file and write header */
1313: vtk_fp = fopen(vtk_filename,"w");
1314: if (!vtk_fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"Cannot open file = %s \n",vtk_filename);
1316: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");
1318: /* coords */
1319: DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1320: N = nx * ny * nz;
1322: #if defined(PETSC_WORDS_BIGENDIAN)
1323: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1324: #else
1325: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1326: #endif
1327: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",si,si+nx-1,sj,sj+ny-1,sk,sk+nz-1);
1328: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",si,si+nx-1,sj,sj+ny-1,sk,sk+nz-1);
1330: memory_offset = 0;
1332: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <CellData></CellData>\n");
1334: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <Points>\n");
1336: /* copy coordinates */
1337: DMGetCoordinateDM(da,&cda);
1338: DMGetCoordinatesLocal(da,&coords);
1339: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%d\" />\n",memory_offset);
1340: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N*3;
1342: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </Points>\n");
1344: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PointData Scalars=\" ");
1345: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_fields,0,0,0,0,0);
1346: for (f=0; f<n_fields; f++) {
1347: const char *field_name;
1348: DMDAGetFieldName(da,f,&field_name);
1349: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"%s ",field_name);
1350: }
1351: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\">\n");
1353: for (f=0; f<n_fields; f++) {
1354: const char *field_name;
1356: DMDAGetFieldName(da,f,&field_name);
1357: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <DataArray type=\"Float64\" Name=\"%s\" format=\"appended\" offset=\"%d\"/>\n", field_name,memory_offset);
1358: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N;
1359: }
1361: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PointData>\n");
1362: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </Piece>\n");
1363: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </StructuredGrid>\n");
1365: PetscMalloc1(N,&buffer);
1366: DMGetLocalVector(da,&l_FIELD);
1367: DMGlobalToLocalBegin(da, FIELD,INSERT_VALUES,l_FIELD);
1368: DMGlobalToLocalEnd(da,FIELD,INSERT_VALUES,l_FIELD);
1369: VecGetArray(l_FIELD,&_L_FIELD);
1371: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <AppendedData encoding=\"raw\">\n");
1372: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"_");
1374: /* write coordinates */
1375: {
1376: int length = sizeof(PetscScalar)*N*3;
1377: PetscScalar *allcoords;
1379: fwrite(&length,sizeof(int),1,vtk_fp);
1380: VecGetArray(coords,&allcoords);
1381: fwrite(allcoords,sizeof(PetscScalar),3*N,vtk_fp);
1382: VecRestoreArray(coords,&allcoords);
1383: }
1384: /* write fields */
1385: for (f=0; f<n_fields; f++) {
1386: int length = sizeof(PetscScalar)*N;
1387: fwrite(&length,sizeof(int),1,vtk_fp);
1388: /* load */
1389: for (i=0; i<N; i++) buffer[i] = _L_FIELD[n_fields*i + f];
1391: /* write */
1392: fwrite(buffer,sizeof(PetscScalar),N,vtk_fp);
1393: }
1394: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\n </AppendedData>\n");
1396: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");
1398: PetscFree(buffer);
1399: VecRestoreArray(l_FIELD,&_L_FIELD);
1400: DMRestoreLocalVector(da,&l_FIELD);
1402: if (vtk_fp) {
1403: fclose(vtk_fp);
1404: vtk_fp = NULL;
1405: }
1407: return(0);
1408: }
1410: PetscErrorCode DAViewVTK_write_PieceExtend(FILE *vtk_fp,PetscInt indent_level,DM da,const char local_file_prefix[])
1411: {
1412: PetscMPIInt size,rank;
1413: MPI_Comm comm;
1414: const PetscInt *lx,*ly,*lz;
1415: PetscInt M,N,P,pM,pN,pP,sum,*olx,*oly,*olz;
1416: PetscInt *osx,*osy,*osz,*oex,*oey,*oez;
1417: PetscInt i,j,k,II,stencil;
1421: /* create file name */
1422: PetscObjectGetComm((PetscObject)da,&comm);
1423: MPI_Comm_size(comm,&size);
1424: MPI_Comm_rank(comm,&rank);
1426: DMDAGetInfo(da,0,&M,&N,&P,&pM,&pN,&pP,0,&stencil,0,0,0,0);
1427: DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
1429: /* generate start,end list */
1430: PetscMalloc1(pM+1,&olx);
1431: PetscMalloc1(pN+1,&oly);
1432: PetscMalloc1(pP+1,&olz);
1433: sum = 0;
1434: for (i=0; i<pM; i++) {
1435: olx[i] = sum;
1436: sum = sum + lx[i];
1437: }
1438: olx[pM] = sum;
1439: sum = 0;
1440: for (i=0; i<pN; i++) {
1441: oly[i] = sum;
1442: sum = sum + ly[i];
1443: }
1444: oly[pN] = sum;
1445: sum = 0;
1446: for (i=0; i<pP; i++) {
1447: olz[i] = sum;
1448: sum = sum + lz[i];
1449: }
1450: olz[pP] = sum;
1452: PetscMalloc1(pM,&osx);
1453: PetscMalloc1(pN,&osy);
1454: PetscMalloc1(pP,&osz);
1455: PetscMalloc1(pM,&oex);
1456: PetscMalloc1(pN,&oey);
1457: PetscMalloc1(pP,&oez);
1458: for (i=0; i<pM; i++) {
1459: osx[i] = olx[i] - stencil;
1460: oex[i] = olx[i] + lx[i] + stencil;
1461: if (osx[i]<0) osx[i]=0;
1462: if (oex[i]>M) oex[i]=M;
1463: }
1465: for (i=0; i<pN; i++) {
1466: osy[i] = oly[i] - stencil;
1467: oey[i] = oly[i] + ly[i] + stencil;
1468: if (osy[i]<0)osy[i]=0;
1469: if (oey[i]>M)oey[i]=N;
1470: }
1471: for (i=0; i<pP; i++) {
1472: osz[i] = olz[i] - stencil;
1473: oez[i] = olz[i] + lz[i] + stencil;
1474: if (osz[i]<0) osz[i]=0;
1475: if (oez[i]>P) oez[i]=P;
1476: }
1478: for (k=0; k<pP; k++) {
1479: for (j=0; j<pN; j++) {
1480: for (i=0; i<pM; i++) {
1481: char name[PETSC_MAX_PATH_LEN];
1482: PetscInt procid = i + j*pM + k*pM*pN; /* convert proc(i,j,k) to pid */
1483: PetscSNPrintf(name,sizeof(name),"subdomain-%s-p%1.4d.vts",local_file_prefix,procid);
1484: for (II=0; II<indent_level; II++) PetscFPrintf(PETSC_COMM_SELF,vtk_fp," ");
1486: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<Piece Extent=\"%d %d %d %d %d %d\" Source=\"%s\"/>\n",
1487: osx[i],oex[i]-1,
1488: osy[j],oey[j]-1,
1489: osz[k],oez[k]-1,name);
1490: }
1491: }
1492: }
1493: PetscFree(olx);
1494: PetscFree(oly);
1495: PetscFree(olz);
1496: PetscFree(osx);
1497: PetscFree(osy);
1498: PetscFree(osz);
1499: PetscFree(oex);
1500: PetscFree(oey);
1501: PetscFree(oez);
1502: return(0);
1503: }
1505: PetscErrorCode DAView_3DVTK_PStructuredGrid(DM da,const char file_prefix[],const char local_file_prefix[])
1506: {
1507: MPI_Comm comm;
1508: PetscMPIInt size,rank;
1509: char vtk_filename[PETSC_MAX_PATH_LEN];
1510: FILE *vtk_fp = NULL;
1511: PetscInt M,N,P,si,sj,sk,nx,ny,nz;
1512: PetscInt i,dofs;
1516: /* only master generates this file */
1517: PetscObjectGetComm((PetscObject)da,&comm);
1518: MPI_Comm_size(comm,&size);
1519: MPI_Comm_rank(comm,&rank);
1521: if (rank != 0) return(0);
1523: /* create file name */
1524: PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"%s.pvts",file_prefix);
1525: vtk_fp = fopen(vtk_filename,"w");
1526: if (!vtk_fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"Cannot open file = %s \n",vtk_filename);
1528: /* (VTK) generate pvts header */
1529: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");
1531: #if defined(PETSC_WORDS_BIGENDIAN)
1532: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1533: #else
1534: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1535: #endif
1537: /* define size of the nodal mesh based on the cell DM */
1538: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,&dofs,0,0,0,0,0);
1539: DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1540: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PStructuredGrid GhostLevel=\"1\" WholeExtent=\"%d %d %d %d %d %d\">\n",0,M-1,0,N-1,0,P-1); /* note overlap = 1 for Q1 */
1542: /* DUMP THE CELL REFERENCES */
1543: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PCellData>\n");
1544: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PCellData>\n");
1546: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PPoints>\n");
1547: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PDataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"%d\"/>\n",NSD);
1548: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PPoints>\n");
1550: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PPointData>\n");
1551: for (i=0; i<dofs; i++) {
1552: const char *fieldname;
1553: DMDAGetFieldName(da,i,&fieldname);
1554: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n",fieldname);
1555: }
1556: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PPointData>\n");
1558: /* write out the parallel information */
1559: DAViewVTK_write_PieceExtend(vtk_fp,2,da,local_file_prefix);
1561: /* close the file */
1562: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PStructuredGrid>\n");
1563: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");
1565: if (vtk_fp) {
1566: fclose(vtk_fp);
1567: vtk_fp = NULL;
1568: }
1569: return(0);
1570: }
1572: PetscErrorCode DAView3DPVTS(DM da, Vec x,const char NAME[])
1573: {
1574: char vts_filename[PETSC_MAX_PATH_LEN];
1575: char pvts_filename[PETSC_MAX_PATH_LEN];
1579: PetscSNPrintf(vts_filename,sizeof(vts_filename),"%s-mesh",NAME);
1580: DAView_3DVTK_StructuredGrid_appended(da,x,vts_filename);
1582: PetscSNPrintf(pvts_filename,sizeof(pvts_filename),"%s-mesh",NAME);
1583: DAView_3DVTK_PStructuredGrid(da,pvts_filename,vts_filename);
1584: return(0);
1585: }
1587: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
1588: {
1590: PetscReal norms[4];
1591: Vec Br,v,w;
1592: Mat A;
1595: KSPGetOperators(ksp,&A,NULL);
1596: MatCreateVecs(A,&w,&v);
1598: KSPBuildResidual(ksp,v,w,&Br);
1600: VecStrideNorm(Br,0,NORM_2,&norms[0]);
1601: VecStrideNorm(Br,1,NORM_2,&norms[1]);
1602: VecStrideNorm(Br,2,NORM_2,&norms[2]);
1603: VecStrideNorm(Br,3,NORM_2,&norms[3]);
1605: VecDestroy(&v);
1606: VecDestroy(&w);
1608: PetscPrintf(PETSC_COMM_WORLD,"%3D KSP Component U,V,W,P residual norm [ %1.12e, %1.12e, %1.12e, %1.12e ]\n",n,norms[0],norms[1],norms[2],norms[3]);
1609: return(0);
1610: }
1612: static PetscErrorCode PCMGSetupViaCoarsen(PC pc,DM da_fine)
1613: {
1614: PetscInt nlevels,k;
1615: PETSC_UNUSED PetscInt finest;
1616: DM *da_list,*daclist;
1617: Mat R;
1618: PetscErrorCode ierr;
1621: nlevels = 1;
1622: PetscOptionsGetInt(NULL,NULL,"-levels",&nlevels,0);
1624: PetscMalloc1(nlevels,&da_list);
1625: for (k=0; k<nlevels; k++) da_list[k] = NULL;
1626: PetscMalloc1(nlevels,&daclist);
1627: for (k=0; k<nlevels; k++) daclist[k] = NULL;
1629: /* finest grid is nlevels - 1 */
1630: finest = nlevels - 1;
1631: daclist[0] = da_fine;
1632: PetscObjectReference((PetscObject)da_fine);
1633: DMCoarsenHierarchy(da_fine,nlevels-1,&daclist[1]);
1634: for (k=0; k<nlevels; k++) {
1635: da_list[k] = daclist[nlevels-1-k];
1636: DMDASetUniformCoordinates(da_list[k],0.0,1.0,0.0,1.0,0.0,1.0);
1637: }
1639: PCMGSetLevels(pc,nlevels,NULL);
1640: PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
1641: PCMGSetGalerkin(pc,PC_MG_GALERKIN_PMAT);
1643: for (k=1; k<nlevels; k++) {
1644: DMCreateInterpolation(da_list[k-1],da_list[k],&R,NULL);
1645: PCMGSetInterpolation(pc,k,R);
1646: MatDestroy(&R);
1647: }
1649: /* tidy up */
1650: for (k=0; k<nlevels; k++) {
1651: DMDestroy(&da_list[k]);
1652: }
1653: PetscFree(da_list);
1654: PetscFree(daclist);
1655: return(0);
1656: }
1658: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx,PetscInt my,PetscInt mz)
1659: {
1660: DM da_Stokes;
1661: PetscInt u_dof,p_dof,dof,stencil_width;
1662: Mat A,B;
1663: DM vel_cda;
1664: Vec vel_coords;
1665: PetscInt p;
1666: Vec f,X;
1667: DMDACoor3d ***_vel_coords;
1668: PetscInt its;
1669: KSP ksp_S;
1670: PetscInt model_definition = 0;
1671: PetscInt ei,ej,ek,sex,sey,sez,Imx,Jmy,Kmz;
1672: CellProperties cell_properties;
1673: PetscBool write_output = PETSC_FALSE;
1677: /* Generate the da for velocity and pressure */
1678: /* Num nodes in each direction is mx+1, my+1, mz+1 */
1679: u_dof = U_DOFS; /* Vx, Vy - velocities */
1680: p_dof = P_DOFS; /* p - pressure */
1681: dof = u_dof+p_dof;
1682: stencil_width = 1;
1683: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1684: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,NULL,&da_Stokes);
1685: DMSetMatType(da_Stokes,MATAIJ);
1686: DMSetFromOptions(da_Stokes);
1687: DMSetUp(da_Stokes);
1688: DMDASetFieldName(da_Stokes,0,"Vx");
1689: DMDASetFieldName(da_Stokes,1,"Vy");
1690: DMDASetFieldName(da_Stokes,2,"Vz");
1691: DMDASetFieldName(da_Stokes,3,"P");
1693: /* unit box [0,1] x [0,1] x [0,1] */
1694: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.0,1.0);
1696: /* create quadrature point info for PDE coefficients */
1697: CellPropertiesCreate(da_Stokes,&cell_properties);
1699: /* interpolate the coordinates to quadrature points */
1700: DMGetCoordinateDM(da_Stokes,&vel_cda);
1701: DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1702: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1703: DMDAGetElementsCorners(da_Stokes,&sex,&sey,&sez);
1704: DMDAGetElementsSizes(da_Stokes,&Imx,&Jmy,&Kmz);
1705: for (ek = sez; ek < sez+Kmz; ek++) {
1706: for (ej = sey; ej < sey+Jmy; ej++) {
1707: for (ei = sex; ei < sex+Imx; ei++) {
1708: /* get coords for the element */
1709: PetscInt ngp;
1710: PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
1711: PetscScalar el_coords[NSD*NODES_PER_EL];
1712: GaussPointCoefficients *cell;
1714: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1715: GetElementCoords3D(_vel_coords,ei,ej,ek,el_coords);
1716: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1718: for (p = 0; p < GAUSS_POINTS; p++) {
1719: PetscScalar xi_p[NSD],Ni_p[NODES_PER_EL];
1720: PetscScalar gp_x,gp_y,gp_z;
1721: PetscInt n;
1723: xi_p[0] = gp_xi[p][0];
1724: xi_p[1] = gp_xi[p][1];
1725: xi_p[2] = gp_xi[p][2];
1726: ShapeFunctionQ13D_Evaluate(xi_p,Ni_p);
1728: gp_x = gp_y = gp_z = 0.0;
1729: for (n = 0; n < NODES_PER_EL; n++) {
1730: gp_x = gp_x+Ni_p[n]*el_coords[NSD*n];
1731: gp_y = gp_y+Ni_p[n]*el_coords[NSD*n+1];
1732: gp_z = gp_z+Ni_p[n]*el_coords[NSD*n+2];
1733: }
1734: cell->gp_coords[NSD*p] = gp_x;
1735: cell->gp_coords[NSD*p+1] = gp_y;
1736: cell->gp_coords[NSD*p+2] = gp_z;
1737: }
1738: }
1739: }
1740: }
1742: PetscOptionsGetInt(NULL,NULL,"-model",&model_definition,NULL);
1744: switch (model_definition) {
1745: case 0: /* isoviscous */
1746: for (ek = sez; ek < sez+Kmz; ek++) {
1747: for (ej = sey; ej < sey+Jmy; ej++) {
1748: for (ei = sex; ei < sex+Imx; ei++) {
1749: GaussPointCoefficients *cell;
1751: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1752: for (p = 0; p < GAUSS_POINTS; p++) {
1753: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1754: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1755: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1757: cell->eta[p] = 1.0;
1759: cell->fx[p] = 0.0*coord_x;
1760: cell->fy[p] = -PetscSinReal(2.2*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1761: cell->fz[p] = 0.0*coord_z;
1762: cell->hc[p] = 0.0;
1763: }
1764: }
1765: }
1766: }
1767: break;
1769: case 1: /* manufactured */
1770: for (ek = sez; ek < sez+Kmz; ek++) {
1771: for (ej = sey; ej < sey+Jmy; ej++) {
1772: for (ei = sex; ei < sex+Imx; ei++) {
1773: PetscReal eta,Fm[NSD],Fc,pos[NSD];
1774: GaussPointCoefficients *cell;
1776: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1777: for (p = 0; p < GAUSS_POINTS; p++) {
1778: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1779: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1780: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1782: pos[0] = coord_x;
1783: pos[1] = coord_y;
1784: pos[2] = coord_z;
1786: evaluate_MS_FrankKamentski(pos,NULL,NULL,&eta,Fm,&Fc);
1787: cell->eta[p] = eta;
1789: cell->fx[p] = Fm[0];
1790: cell->fy[p] = Fm[1];
1791: cell->fz[p] = Fm[2];
1792: cell->hc[p] = Fc;
1793: }
1794: }
1795: }
1796: }
1797: break;
1799: case 2: /* solcx */
1800: for (ek = sez; ek < sez+Kmz; ek++) {
1801: for (ej = sey; ej < sey+Jmy; ej++) {
1802: for (ei = sex; ei < sex+Imx; ei++) {
1803: GaussPointCoefficients *cell;
1805: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1806: for (p = 0; p < GAUSS_POINTS; p++) {
1807: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1808: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1809: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1811: cell->eta[p] = 1.0;
1813: cell->fx[p] = 0.0;
1814: cell->fy[p] = -PetscSinReal(3.0*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1815: cell->fz[p] = 0.0*coord_z;
1816: cell->hc[p] = 0.0;
1817: }
1818: }
1819: }
1820: }
1821: break;
1823: case 3: /* sinker */
1824: for (ek = sez; ek < sez+Kmz; ek++) {
1825: for (ej = sey; ej < sey+Jmy; ej++) {
1826: for (ei = sex; ei < sex+Imx; ei++) {
1827: GaussPointCoefficients *cell;
1829: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1830: for (p = 0; p < GAUSS_POINTS; p++) {
1831: PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p]);
1832: PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);
1833: PetscReal zp = PetscRealPart(cell->gp_coords[NSD*p+2]);
1835: cell->eta[p] = 1.0e-2;
1836: cell->fx[p] = 0.0;
1837: cell->fy[p] = 0.0;
1838: cell->fz[p] = 0.0;
1839: cell->hc[p] = 0.0;
1841: if ((PetscAbs(xp-0.5) < 0.2) && (PetscAbs(yp-0.5) < 0.2) && (PetscAbs(zp-0.5) < 0.2)) {
1842: cell->eta[p] = 1.0;
1843: cell->fz[p] = 1.0;
1844: }
1846: }
1847: }
1848: }
1849: }
1850: break;
1852: case 4: /* subdomain jumps */
1853: for (ek = sez; ek < sez+Kmz; ek++) {
1854: for (ej = sey; ej < sey+Jmy; ej++) {
1855: for (ei = sex; ei < sex+Imx; ei++) {
1856: PetscReal opts_mag,opts_eta0;
1857: PetscInt px,py,pz;
1858: PetscBool jump;
1859: PetscMPIInt rr;
1860: GaussPointCoefficients *cell;
1862: opts_mag = 1.0;
1863: opts_eta0 = 1.e-2;
1865: PetscOptionsGetReal(NULL,NULL,"-jump_eta0",&opts_eta0,NULL);
1866: PetscOptionsGetReal(NULL,NULL,"-jump_magnitude",&opts_mag,NULL);
1867: DMDAGetInfo(da_Stokes,NULL,NULL,NULL,NULL,&px,&py,&pz,NULL,NULL,NULL,NULL,NULL,NULL);
1868: rr = PetscGlobalRank%(px*py);
1869: if (px%2) jump = (PetscBool)(rr%2);
1870: else jump = (PetscBool)((rr/px)%2 ? rr%2 : !(rr%2));
1871: rr = PetscGlobalRank/(px*py);
1872: if (rr%2) jump = (PetscBool)!jump;
1873: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1874: for (p = 0; p < GAUSS_POINTS; p++) {
1875: PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p]);
1876: PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);
1878: cell->eta[p] = jump ? PetscPowReal(10.0,opts_mag) : opts_eta0;
1879: cell->fx[p] = 0.0;
1880: cell->fy[p] = -PetscSinReal(2.2*PETSC_PI*yp)*PetscCosReal(1.0*PETSC_PI*xp);
1881: cell->fz[p] = 0.0;
1882: cell->hc[p] = 0.0;
1883: }
1884: }
1885: }
1886: }
1887: break;
1888: default:
1889: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"No default model is supported. Choose either -model {0,1,2,3}");
1890: break;
1891: }
1893: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1895: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1896: DMCreateMatrix(da_Stokes,&A);
1897: DMCreateMatrix(da_Stokes,&B);
1898: DMCreateGlobalVector(da_Stokes,&X);
1899: DMCreateGlobalVector(da_Stokes,&f);
1901: /* assemble A11 */
1902: MatZeroEntries(A);
1903: MatZeroEntries(B);
1904: VecZeroEntries(f);
1906: AssembleA_Stokes(A,da_Stokes,cell_properties);
1907: MatViewFromOptions(A,NULL,"-amat_view");
1908: AssembleA_PCStokes(B,da_Stokes,cell_properties);
1909: MatViewFromOptions(B,NULL,"-bmat_view");
1910: /* build force vector */
1911: AssembleF_Stokes(f,da_Stokes,cell_properties);
1913: /* SOLVE */
1914: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1915: KSPSetOptionsPrefix(ksp_S,"stokes_"); /* stokes */
1916: KSPSetOperators(ksp_S,A,B);
1917: KSPSetFromOptions(ksp_S);
1919: {
1920: PC pc;
1921: const PetscInt ufields[] = {0,1,2},pfields[] = {3};
1922: KSPGetPC(ksp_S,&pc);
1923: PCFieldSplitSetBlockSize(pc,4);
1924: PCFieldSplitSetFields(pc,"u",3,ufields,ufields);
1925: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1926: }
1928: {
1929: PC pc;
1930: PetscBool same = PETSC_FALSE;
1931: KSPGetPC(ksp_S,&pc);
1932: PetscObjectTypeCompare((PetscObject)pc,PCMG,&same);
1933: if (same) {
1934: PCMGSetupViaCoarsen(pc,da_Stokes);
1935: }
1936: }
1938: {
1939: PC pc;
1940: PetscBool same = PETSC_FALSE;
1941: KSPGetPC(ksp_S,&pc);
1942: PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&same);
1943: if (same) {
1944: KSPSetOperators(ksp_S,A,A);
1945: }
1946: }
1948: {
1949: PetscBool stokes_monitor = PETSC_FALSE;
1950: PetscOptionsGetBool(NULL,NULL,"-stokes_ksp_monitor_blocks",&stokes_monitor,0);
1951: if (stokes_monitor) {
1952: KSPMonitorSet(ksp_S,KSPMonitorStokesBlocks,NULL,NULL);
1953: }
1954: }
1955: KSPSolve(ksp_S,f,X);
1957: PetscOptionsGetBool(NULL,NULL,"-write_pvts",&write_output,NULL);
1958: if (write_output) {DAView3DPVTS(da_Stokes,X,"up");}
1959: {
1960: PetscBool flg = PETSC_FALSE;
1961: char filename[PETSC_MAX_PATH_LEN];
1962: PetscOptionsGetString(NULL,NULL,"-write_binary",filename,sizeof(filename),&flg);
1963: if (flg) {
1964: PetscViewer viewer;
1965: /* PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer); */
1966: PetscViewerVTKOpen(PETSC_COMM_WORLD,"ex42.vts",FILE_MODE_WRITE,&viewer);
1967: VecView(X,viewer);
1968: PetscViewerDestroy(&viewer);
1969: }
1970: }
1971: KSPGetIterationNumber(ksp_S,&its);
1973: /* verify */
1974: if (model_definition == 1) {
1975: DM da_Stokes_analytic;
1976: Vec X_analytic;
1978: DMDACreateManufacturedSolution(mx,my,mz,&da_Stokes_analytic,&X_analytic);
1979: if (write_output) {
1980: DAView3DPVTS(da_Stokes_analytic,X_analytic,"ms");
1981: }
1982: DMDAIntegrateErrors3D(da_Stokes_analytic,X,X_analytic);
1983: if (write_output) {
1984: DAView3DPVTS(da_Stokes,X,"up2");
1985: }
1986: DMDestroy(&da_Stokes_analytic);
1987: VecDestroy(&X_analytic);
1988: }
1990: KSPDestroy(&ksp_S);
1991: VecDestroy(&X);
1992: VecDestroy(&f);
1993: MatDestroy(&A);
1994: MatDestroy(&B);
1996: CellPropertiesDestroy(&cell_properties);
1997: DMDestroy(&da_Stokes);
1998: return(0);
1999: }
2001: int main(int argc,char **args)
2002: {
2004: PetscInt mx,my,mz;
2006: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
2007: mx = my = mz = 10;
2008: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
2009: my = mx; mz = mx;
2010: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
2011: PetscOptionsGetInt(NULL,NULL,"-mz",&mz,NULL);
2012: solve_stokes_3d_coupled(mx,my,mz);
2013: PetscFinalize();
2014: return ierr;
2015: }
2018: /*TEST
2020: test:
2021: args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type lu
2023: test:
2024: suffix: 2
2025: nsize: 3
2026: args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type redundant -stokes_redundant_pc_type lu
2028: test:
2029: suffix: bddc_stokes
2030: nsize: 6
2031: args: -mx 5 -my 4 -mz 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_dirichlet_pc_type svd -stokes_pc_bddc_neumann_pc_type svd -stokes_pc_bddc_coarse_redundant_pc_type svd
2033: test:
2034: suffix: bddc_stokes_deluxe
2035: nsize: 6
2036: args: -mx 5 -my 4 -mz 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_dirichlet_pc_type svd -stokes_pc_bddc_neumann_pc_type svd -stokes_pc_bddc_coarse_redundant_pc_type svd -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc
2038: test:
2039: requires: !single
2040: suffix: bddc_stokes_subdomainjump_deluxe
2041: nsize: 9
2042: args: -model 4 -jump_magnitude 4 -mx 6 -my 6 -mz 2 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc -stokes_pc_bddc_schur_layers 1
2045: TEST*/