Actual source code: ex42.c
petsc-3.8.4 2018-03-24
1: static char help[] = "Solves the incompressible, variable viscosity stokes equation in 3d using Q1Q1 elements, \n\
2: stabilized with Bochev's polynomial projection method. Note that implementation here assumes \n\
3: all boundaries are free-slip, i.e. zero normal flow and zero tangential stress \n\
4: -mx : number elements in x-direction \n\
5: -my : number elements in y-direction \n\
6: -mz : number elements in z-direction \n\
7: -stokes_ksp_monitor_blocks : active monitor for each component u,v,w,p \n\
8: -model : defines viscosity and forcing function. Choose either: 0 (isoviscous), 1 (manufactured-broken), 2 (solcx-2d), 3 (sinker) \n";
10: /* Contributed by Dave May */
12: #include <petscksp.h>
13: #include <petscdmda.h>
15: #define PROFILE_TIMING
16: #define ASSEMBLE_LOWER_TRIANGULAR
18: #define NSD 3 /* number of spatial dimensions */
19: #define NODES_PER_EL 8 /* nodes per element */
20: #define U_DOFS 3 /* degrees of freedom per velocity node */
21: #define P_DOFS 1 /* degrees of freedom per pressure node */
22: #define GAUSS_POINTS 8
24: /* Gauss point based evaluation */
25: typedef struct {
26: PetscScalar gp_coords[NSD*GAUSS_POINTS];
27: PetscScalar eta[GAUSS_POINTS];
28: PetscScalar fx[GAUSS_POINTS];
29: PetscScalar fy[GAUSS_POINTS];
30: PetscScalar fz[GAUSS_POINTS];
31: PetscScalar hc[GAUSS_POINTS];
32: } GaussPointCoefficients;
34: typedef struct {
35: PetscScalar u_dof;
36: PetscScalar v_dof;
37: PetscScalar w_dof;
38: PetscScalar p_dof;
39: } StokesDOF;
41: typedef struct _p_CellProperties *CellProperties;
42: struct _p_CellProperties {
43: PetscInt ncells;
44: PetscInt mx,my,mz;
45: PetscInt sex,sey,sez;
46: GaussPointCoefficients *gpc;
47: };
49: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz);
51: /* elements */
52: PetscErrorCode CellPropertiesCreate(DM da_stokes,CellProperties *C)
53: {
55: CellProperties cells;
56: PetscInt mx,my,mz,sex,sey,sez;
59: PetscNew(&cells);
61: DMDAGetElementCorners(da_stokes,&sex,&sey,&sez,&mx,&my,&mz);
63: cells->mx = mx;
64: cells->my = my;
65: cells->mz = mz;
66: cells->ncells = mx * my * mz;
67: cells->sex = sex;
68: cells->sey = sey;
69: cells->sez = sez;
71: PetscMalloc1(mx*my*mz,&cells->gpc);
73: *C = cells;
74: return(0);
75: }
77: PetscErrorCode CellPropertiesDestroy(CellProperties *C)
78: {
80: CellProperties cells;
83: if (!C) return(0);
84: cells = *C;
85: PetscFree(cells->gpc);
86: PetscFree(cells);
87: *C = NULL;
88: return(0);
89: }
91: PetscErrorCode CellPropertiesGetCell(CellProperties C,PetscInt II,PetscInt J,PetscInt K,GaussPointCoefficients **G)
92: {
94: *G = &C->gpc[(II-C->sex) + (J-C->sey)*C->mx + (K-C->sez)*C->mx*C->my];
95: return(0);
96: }
98: /* FEM routines */
99: /*
100: Element: Local basis function ordering
101: 1-----2
102: | |
103: | |
104: 0-----3
105: */
106: static void ShapeFunctionQ13D_Evaluate(PetscScalar _xi[],PetscScalar Ni[])
107: {
108: PetscReal xi = PetscRealPart(_xi[0]);
109: PetscReal eta = PetscRealPart(_xi[1]);
110: PetscReal zeta = PetscRealPart(_xi[2]);
112: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
113: Ni[1] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
114: Ni[2] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
115: Ni[3] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
117: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
118: Ni[5] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
119: Ni[6] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
120: Ni[7] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
121: }
123: static void ShapeFunctionQ13D_Evaluate_dxi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
124: {
125: PetscReal xi = PetscRealPart(_xi[0]);
126: PetscReal eta = PetscRealPart(_xi[1]);
127: PetscReal zeta = PetscRealPart(_xi[2]);
128: /* xi deriv */
129: GNi[0][0] = -0.125 * (1.0 - eta) * (1.0 - zeta);
130: GNi[0][1] = -0.125 * (1.0 + eta) * (1.0 - zeta);
131: GNi[0][2] = 0.125 * (1.0 + eta) * (1.0 - zeta);
132: GNi[0][3] = 0.125 * (1.0 - eta) * (1.0 - zeta);
134: GNi[0][4] = -0.125 * (1.0 - eta) * (1.0 + zeta);
135: GNi[0][5] = -0.125 * (1.0 + eta) * (1.0 + zeta);
136: GNi[0][6] = 0.125 * (1.0 + eta) * (1.0 + zeta);
137: GNi[0][7] = 0.125 * (1.0 - eta) * (1.0 + zeta);
138: /* eta deriv */
139: GNi[1][0] = -0.125 * (1.0 - xi) * (1.0 - zeta);
140: GNi[1][1] = 0.125 * (1.0 - xi) * (1.0 - zeta);
141: GNi[1][2] = 0.125 * (1.0 + xi) * (1.0 - zeta);
142: GNi[1][3] = -0.125 * (1.0 + xi) * (1.0 - zeta);
144: GNi[1][4] = -0.125 * (1.0 - xi) * (1.0 + zeta);
145: GNi[1][5] = 0.125 * (1.0 - xi) * (1.0 + zeta);
146: GNi[1][6] = 0.125 * (1.0 + xi) * (1.0 + zeta);
147: GNi[1][7] = -0.125 * (1.0 + xi) * (1.0 + zeta);
148: /* zeta deriv */
149: GNi[2][0] = -0.125 * (1.0 - xi) * (1.0 - eta);
150: GNi[2][1] = -0.125 * (1.0 - xi) * (1.0 + eta);
151: GNi[2][2] = -0.125 * (1.0 + xi) * (1.0 + eta);
152: GNi[2][3] = -0.125 * (1.0 + xi) * (1.0 - eta);
154: GNi[2][4] = 0.125 * (1.0 - xi) * (1.0 - eta);
155: GNi[2][5] = 0.125 * (1.0 - xi) * (1.0 + eta);
156: GNi[2][6] = 0.125 * (1.0 + xi) * (1.0 + eta);
157: GNi[2][7] = 0.125 * (1.0 + xi) * (1.0 - eta);
158: }
160: static void matrix_inverse_3x3(PetscScalar A[3][3],PetscScalar B[3][3])
161: {
162: PetscScalar t4, t6, t8, t10, t12, t14, t17;
164: t4 = A[2][0] * A[0][1];
165: t6 = A[2][0] * A[0][2];
166: t8 = A[1][0] * A[0][1];
167: t10 = A[1][0] * A[0][2];
168: t12 = A[0][0] * A[1][1];
169: t14 = A[0][0] * A[1][2];
170: t17 = 0.1e1 / (t4 * A[1][2] - t6 * A[1][1] - t8 * A[2][2] + t10 * A[2][1] + t12 * A[2][2] - t14 * A[2][1]);
172: B[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) * t17;
173: B[0][1] = -(A[0][1] * A[2][2] - A[0][2] * A[2][1]) * t17;
174: B[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) * t17;
175: B[1][0] = -(-A[2][0] * A[1][2] + A[1][0] * A[2][2]) * t17;
176: B[1][1] = (-t6 + A[0][0] * A[2][2]) * t17;
177: B[1][2] = -(-t10 + t14) * t17;
178: B[2][0] = (-A[2][0] * A[1][1] + A[1][0] * A[2][1]) * t17;
179: B[2][1] = -(-t4 + A[0][0] * A[2][1]) * t17;
180: B[2][2] = (-t8 + t12) * t17;
181: }
183: static void ShapeFunctionQ13D_Evaluate_dx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
184: {
185: PetscScalar J00,J01,J02,J10,J11,J12,J20,J21,J22;
186: PetscInt n;
187: PetscScalar iJ[3][3],JJ[3][3];
189: J00 = J01 = J02 = 0.0;
190: J10 = J11 = J12 = 0.0;
191: J20 = J21 = J22 = 0.0;
192: for (n=0; n<NODES_PER_EL; n++) {
193: PetscScalar cx = coords[NSD*n + 0];
194: PetscScalar cy = coords[NSD*n + 1];
195: PetscScalar cz = coords[NSD*n + 2];
197: /* J_ij = d(x_j) / d(xi_i) */ /* J_ij = \sum _I GNi[j][I} * x_i */
198: J00 = J00 + GNi[0][n] * cx; /* J_xx */
199: J01 = J01 + GNi[0][n] * cy; /* J_xy = dx/deta */
200: J02 = J02 + GNi[0][n] * cz; /* J_xz = dx/dzeta */
202: J10 = J10 + GNi[1][n] * cx; /* J_yx = dy/dxi */
203: J11 = J11 + GNi[1][n] * cy; /* J_yy */
204: J12 = J12 + GNi[1][n] * cz; /* J_yz */
206: J20 = J20 + GNi[2][n] * cx; /* J_zx */
207: J21 = J21 + GNi[2][n] * cy; /* J_zy */
208: J22 = J22 + GNi[2][n] * cz; /* J_zz */
209: }
211: JJ[0][0] = J00; JJ[0][1] = J01; JJ[0][2] = J02;
212: JJ[1][0] = J10; JJ[1][1] = J11; JJ[1][2] = J12;
213: JJ[2][0] = J20; JJ[2][1] = J21; JJ[2][2] = J22;
215: matrix_inverse_3x3(JJ,iJ);
217: *det_J = J00*J11*J22 - J00*J12*J21 - J10*J01*J22 + J10*J02*J21 + J20*J01*J12 - J20*J02*J11;
219: for (n=0; n<NODES_PER_EL; n++) {
220: GNx[0][n] = GNi[0][n]*iJ[0][0] + GNi[1][n]*iJ[0][1] + GNi[2][n]*iJ[0][2];
221: GNx[1][n] = GNi[0][n]*iJ[1][0] + GNi[1][n]*iJ[1][1] + GNi[2][n]*iJ[1][2];
222: GNx[2][n] = GNi[0][n]*iJ[2][0] + GNi[1][n]*iJ[2][1] + GNi[2][n]*iJ[2][2];
223: }
224: }
226: static void ConstructGaussQuadrature3D(PetscInt *ngp,PetscScalar gp_xi[][NSD],PetscScalar gp_weight[])
227: {
228: *ngp = 8;
229: gp_xi[0][0] = -0.57735026919; gp_xi[0][1] = -0.57735026919; gp_xi[0][2] = -0.57735026919;
230: gp_xi[1][0] = -0.57735026919; gp_xi[1][1] = 0.57735026919; gp_xi[1][2] = -0.57735026919;
231: gp_xi[2][0] = 0.57735026919; gp_xi[2][1] = 0.57735026919; gp_xi[2][2] = -0.57735026919;
232: gp_xi[3][0] = 0.57735026919; gp_xi[3][1] = -0.57735026919; gp_xi[3][2] = -0.57735026919;
234: gp_xi[4][0] = -0.57735026919; gp_xi[4][1] = -0.57735026919; gp_xi[4][2] = 0.57735026919;
235: gp_xi[5][0] = -0.57735026919; gp_xi[5][1] = 0.57735026919; gp_xi[5][2] = 0.57735026919;
236: gp_xi[6][0] = 0.57735026919; gp_xi[6][1] = 0.57735026919; gp_xi[6][2] = 0.57735026919;
237: gp_xi[7][0] = 0.57735026919; gp_xi[7][1] = -0.57735026919; gp_xi[7][2] = 0.57735026919;
239: gp_weight[0] = 1.0;
240: gp_weight[1] = 1.0;
241: gp_weight[2] = 1.0;
242: gp_weight[3] = 1.0;
244: gp_weight[4] = 1.0;
245: gp_weight[5] = 1.0;
246: gp_weight[6] = 1.0;
247: gp_weight[7] = 1.0;
248: }
250: /* procs to the left claim the ghost node as their element */
251: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
252: {
253: PetscInt m,n,p,M,N,P;
254: PetscInt sx,sy,sz;
258: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
259: DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);
261: if (mxl) {
262: *mxl = m;
263: if ((sx+m) == M) *mxl = m-1; /* last proc */
264: }
265: if (myl) {
266: *myl = n;
267: if ((sy+n) == N) *myl = n-1; /* last proc */
268: }
269: if (mzl) {
270: *mzl = p;
271: if ((sz+p) == P) *mzl = p-1; /* last proc */
272: }
273: return(0);
274: }
276: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
277: {
278: PetscInt si,sj,sk;
282: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
284: if (sx) {
285: *sx = si;
286: if (si != 0) *sx = si+1;
287: }
288: if (sy) {
289: *sy = sj;
290: if (sj != 0) *sy = sj+1;
291: }
292: if (sz) {
293: *sz = sk;
294: if (sk != 0) *sz = sk+1;
295: }
296: DMDAGetLocalElementSize(da,mx,my,mz);
297: return(0);
298: }
300: /*
301: i,j are the element indices
302: The unknown is a vector quantity.
303: The s[].c is used to indicate the degree of freedom.
304: */
305: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j,PetscInt k)
306: {
307: PetscInt n;
310: /* velocity */
311: n = 0;
312: /* node 0 */
313: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++; /* Vx0 */
314: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++; /* Vy0 */
315: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++; /* Vz0 */
317: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
318: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
319: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
321: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
322: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
323: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
325: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++;
326: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++;
327: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++;
329: /* */
330: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++; /* Vx4 */
331: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++; /* Vy4 */
332: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++; /* Vz4 */
334: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
335: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
336: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
338: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
339: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
340: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
342: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++;
343: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++;
344: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++;
346: /* pressure */
347: n = 0;
349: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++; /* P0 */
350: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
351: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
352: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++;
354: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++; /* P0 */
355: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
356: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
357: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++;
358: return(0);
359: }
361: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords,PetscInt i,PetscInt j,PetscInt k,PetscScalar el_coord[])
362: {
364: /* get coords for the element */
365: el_coord[0] = coords[k][j][i].x;
366: el_coord[1] = coords[k][j][i].y;
367: el_coord[2] = coords[k][j][i].z;
369: el_coord[3] = coords[k][j+1][i].x;
370: el_coord[4] = coords[k][j+1][i].y;
371: el_coord[5] = coords[k][j+1][i].z;
373: el_coord[6] = coords[k][j+1][i+1].x;
374: el_coord[7] = coords[k][j+1][i+1].y;
375: el_coord[8] = coords[k][j+1][i+1].z;
377: el_coord[9] = coords[k][j][i+1].x;
378: el_coord[10] = coords[k][j][i+1].y;
379: el_coord[11] = coords[k][j][i+1].z;
381: el_coord[12] = coords[k+1][j][i].x;
382: el_coord[13] = coords[k+1][j][i].y;
383: el_coord[14] = coords[k+1][j][i].z;
385: el_coord[15] = coords[k+1][j+1][i].x;
386: el_coord[16] = coords[k+1][j+1][i].y;
387: el_coord[17] = coords[k+1][j+1][i].z;
389: el_coord[18] = coords[k+1][j+1][i+1].x;
390: el_coord[19] = coords[k+1][j+1][i+1].y;
391: el_coord[20] = coords[k+1][j+1][i+1].z;
393: el_coord[21] = coords[k+1][j][i+1].x;
394: el_coord[22] = coords[k+1][j][i+1].y;
395: el_coord[23] = coords[k+1][j][i+1].z;
396: return(0);
397: }
399: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field,PetscInt i,PetscInt j,PetscInt k,StokesDOF nodal_fields[])
400: {
402: /* get the nodal fields for u */
403: nodal_fields[0].u_dof = field[k][j][i].u_dof;
404: nodal_fields[0].v_dof = field[k][j][i].v_dof;
405: nodal_fields[0].w_dof = field[k][j][i].w_dof;
407: nodal_fields[1].u_dof = field[k][j+1][i].u_dof;
408: nodal_fields[1].v_dof = field[k][j+1][i].v_dof;
409: nodal_fields[1].w_dof = field[k][j+1][i].w_dof;
411: nodal_fields[2].u_dof = field[k][j+1][i+1].u_dof;
412: nodal_fields[2].v_dof = field[k][j+1][i+1].v_dof;
413: nodal_fields[2].w_dof = field[k][j+1][i+1].w_dof;
415: nodal_fields[3].u_dof = field[k][j][i+1].u_dof;
416: nodal_fields[3].v_dof = field[k][j][i+1].v_dof;
417: nodal_fields[3].w_dof = field[k][j][i+1].w_dof;
419: nodal_fields[4].u_dof = field[k+1][j][i].u_dof;
420: nodal_fields[4].v_dof = field[k+1][j][i].v_dof;
421: nodal_fields[4].w_dof = field[k+1][j][i].w_dof;
423: nodal_fields[5].u_dof = field[k+1][j+1][i].u_dof;
424: nodal_fields[5].v_dof = field[k+1][j+1][i].v_dof;
425: nodal_fields[5].w_dof = field[k+1][j+1][i].w_dof;
427: nodal_fields[6].u_dof = field[k+1][j+1][i+1].u_dof;
428: nodal_fields[6].v_dof = field[k+1][j+1][i+1].v_dof;
429: nodal_fields[6].w_dof = field[k+1][j+1][i+1].w_dof;
431: nodal_fields[7].u_dof = field[k+1][j][i+1].u_dof;
432: nodal_fields[7].v_dof = field[k+1][j][i+1].v_dof;
433: nodal_fields[7].w_dof = field[k+1][j][i+1].w_dof;
435: /* get the nodal fields for p */
436: nodal_fields[0].p_dof = field[k][j][i].p_dof;
437: nodal_fields[1].p_dof = field[k][j+1][i].p_dof;
438: nodal_fields[2].p_dof = field[k][j+1][i+1].p_dof;
439: nodal_fields[3].p_dof = field[k][j][i+1].p_dof;
441: nodal_fields[4].p_dof = field[k+1][j][i].p_dof;
442: nodal_fields[5].p_dof = field[k+1][j+1][i].p_dof;
443: nodal_fields[6].p_dof = field[k+1][j+1][i+1].p_dof;
444: nodal_fields[7].p_dof = field[k+1][j][i+1].p_dof;
445: return(0);
446: }
448: static PetscInt ASS_MAP_wIwDI_uJuDJ(PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)
449: {
450: PetscInt ij;
451: PETSC_UNUSED PetscInt r,c,nr,nc;
453: nr = w_NPE*w_dof;
454: nc = u_NPE*u_dof;
456: r = w_dof*wi+wd;
457: c = u_dof*ui+ud;
459: ij = r*nc+c;
461: return ij;
462: }
464: static PetscErrorCode DMDASetValuesLocalStencil3D_ADD_VALUES(StokesDOF ***fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
465: {
466: PetscInt n,II,J,K;
469: for (n = 0; n<NODES_PER_EL; n++) {
470: II = u_eqn[NSD*n].i;
471: J = u_eqn[NSD*n].j;
472: K = u_eqn[NSD*n].k;
474: fields_F[K][J][II].u_dof = fields_F[K][J][II].u_dof+Fe_u[NSD*n];
476: II = u_eqn[NSD*n+1].i;
477: J = u_eqn[NSD*n+1].j;
478: K = u_eqn[NSD*n+1].k;
480: fields_F[K][J][II].v_dof = fields_F[K][J][II].v_dof+Fe_u[NSD*n+1];
482: II = u_eqn[NSD*n+2].i;
483: J = u_eqn[NSD*n+2].j;
484: K = u_eqn[NSD*n+2].k;
485: fields_F[K][J][II].w_dof = fields_F[K][J][II].w_dof+Fe_u[NSD*n+2];
487: II = p_eqn[n].i;
488: J = p_eqn[n].j;
489: K = p_eqn[n].k;
491: fields_F[K][J][II].p_dof = fields_F[K][J][II].p_dof+Fe_p[n];
493: }
494: return(0);
495: }
497: static void FormStressOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
498: {
499: PetscInt ngp;
500: PetscScalar gp_xi[GAUSS_POINTS][NSD];
501: PetscScalar gp_weight[GAUSS_POINTS];
502: PetscInt p,i,j,k;
503: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
504: PetscScalar J_p,tildeD[6];
505: PetscScalar B[6][U_DOFS*NODES_PER_EL];
506: const PetscInt nvdof = U_DOFS*NODES_PER_EL;
508: /* define quadrature rule */
509: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
511: /* evaluate integral */
512: for (p = 0; p < ngp; p++) {
513: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
514: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
516: for (i = 0; i < NODES_PER_EL; i++) {
517: PetscScalar d_dx_i = GNx_p[0][i];
518: PetscScalar d_dy_i = GNx_p[1][i];
519: PetscScalar d_dz_i = GNx_p[2][i];
521: B[0][3*i] = d_dx_i; B[0][3*i+1] = 0.0; B[0][3*i+2] = 0.0;
522: B[1][3*i] = 0.0; B[1][3*i+1] = d_dy_i; B[1][3*i+2] = 0.0;
523: B[2][3*i] = 0.0; B[2][3*i+1] = 0.0; B[2][3*i+2] = d_dz_i;
525: B[3][3*i] = d_dy_i; B[3][3*i+1] = d_dx_i; B[3][3*i+2] = 0.0; /* e_xy */
526: B[4][3*i] = d_dz_i; B[4][3*i+1] = 0.0; B[4][3*i+2] = d_dx_i; /* e_xz */
527: B[5][3*i] = 0.0; B[5][3*i+1] = d_dz_i; B[5][3*i+2] = d_dy_i; /* e_yz */
528: }
531: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
532: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
533: tildeD[2] = 2.0*gp_weight[p]*J_p*eta[p];
535: tildeD[3] = gp_weight[p]*J_p*eta[p];
536: tildeD[4] = gp_weight[p]*J_p*eta[p];
537: tildeD[5] = gp_weight[p]*J_p*eta[p];
539: /* form Bt tildeD B */
540: /*
541: Ke_ij = Bt_ik . D_kl . B_lj
542: = B_ki . D_kl . B_lj
543: = B_ki . D_kk . B_kj
544: */
545: for (i = 0; i < nvdof; i++) {
546: for (j = i; j < nvdof; j++) {
547: for (k = 0; k < 6; k++) {
548: Ke[i*nvdof+j] += B[k][i]*tildeD[k]*B[k][j];
549: }
550: }
551: }
553: }
554: /* fill lower triangular part */
555: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
556: for (i = 0; i < nvdof; i++) {
557: for (j = i; j < nvdof; j++) {
558: Ke[j*nvdof+i] = Ke[i*nvdof+j];
559: }
560: }
561: #endif
562: }
564: static void FormGradientOperatorQ13D(PetscScalar Ke[],PetscScalar coords[])
565: {
566: PetscInt ngp;
567: PetscScalar gp_xi[GAUSS_POINTS][NSD];
568: PetscScalar gp_weight[GAUSS_POINTS];
569: PetscInt p,i,j,di;
570: PetscScalar Ni_p[NODES_PER_EL];
571: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
572: PetscScalar J_p,fac;
574: /* define quadrature rule */
575: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
577: /* evaluate integral */
578: for (p = 0; p < ngp; p++) {
579: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
580: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
581: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
582: fac = gp_weight[p]*J_p;
584: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
585: for (di = 0; di < NSD; di++) { /* u dofs */
586: for (j = 0; j < NODES_PER_EL; j++) { /* p nodes, p dofs = 1 (ie no loop) */
587: PetscInt IJ;
588: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,3,j,0,NODES_PER_EL,1);
590: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
591: }
592: }
593: }
594: }
595: }
597: static void FormDivergenceOperatorQ13D(PetscScalar De[],PetscScalar coords[])
598: {
599: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
600: PetscInt i,j;
601: PetscInt nr_g,nc_g;
603: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
604: FormGradientOperatorQ13D(Ge,coords);
606: nr_g = U_DOFS*NODES_PER_EL;
607: nc_g = P_DOFS*NODES_PER_EL;
609: for (i = 0; i < nr_g; i++) {
610: for (j = 0; j < nc_g; j++) {
611: De[nr_g*j+i] = Ge[nc_g*i+j];
612: }
613: }
614: }
616: static void FormStabilisationOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
617: {
618: PetscInt ngp;
619: PetscScalar gp_xi[GAUSS_POINTS][NSD];
620: PetscScalar gp_weight[GAUSS_POINTS];
621: PetscInt p,i,j;
622: PetscScalar Ni_p[NODES_PER_EL];
623: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
624: PetscScalar J_p,fac,eta_avg;
626: /* define quadrature rule */
627: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
629: /* evaluate integral */
630: for (p = 0; p < ngp; p++) {
631: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
632: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
633: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
634: fac = gp_weight[p]*J_p;
635: /*
636: for (i = 0; i < NODES_PER_EL; i++) {
637: for (j = i; j < NODES_PER_EL; j++) {
638: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
639: }
640: }
641: */
643: for (i = 0; i < NODES_PER_EL; i++) {
644: for (j = 0; j < NODES_PER_EL; j++) {
645: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
646: }
647: }
648: }
650: /* scale */
651: eta_avg = 0.0;
652: for (p = 0; p < ngp; p++) eta_avg += eta[p];
653: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
654: fac = 1.0/eta_avg;
656: /*
657: for (i = 0; i < NODES_PER_EL; i++) {
658: for (j = i; j < NODES_PER_EL; j++) {
659: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
660: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
661: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
662: #endif
663: }
664: }
665: */
667: for (i = 0; i < NODES_PER_EL; i++) {
668: for (j = 0; j < NODES_PER_EL; j++) {
669: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
670: }
671: }
672: }
674: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
675: {
676: PetscInt ngp;
677: PetscScalar gp_xi[GAUSS_POINTS][NSD];
678: PetscScalar gp_weight[GAUSS_POINTS];
679: PetscInt p,i,j;
680: PetscScalar Ni_p[NODES_PER_EL];
681: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
682: PetscScalar J_p,fac,eta_avg;
684: /* define quadrature rule */
685: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
687: /* evaluate integral */
688: for (p = 0; p < ngp; p++) {
689: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
690: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
691: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
692: fac = gp_weight[p]*J_p;
694: /*
695: for (i = 0; i < NODES_PER_EL; i++) {
696: for (j = i; j < NODES_PER_EL; j++) {
697: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
698: }
699: }
700: */
702: for (i = 0; i < NODES_PER_EL; i++) {
703: for (j = 0; j < NODES_PER_EL; j++) {
704: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
705: }
706: }
707: }
709: /* scale */
710: eta_avg = 0.0;
711: for (p = 0; p < ngp; p++) eta_avg += eta[p];
712: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
713: fac = 1.0/eta_avg;
714: /*
715: for (i = 0; i < NODES_PER_EL; i++) {
716: for (j = i; j < NODES_PER_EL; j++) {
717: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
718: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
719: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
720: #endif
721: }
722: }
723: */
725: for (i = 0; i < NODES_PER_EL; i++) {
726: for (j = 0; j < NODES_PER_EL; j++) {
727: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
728: }
729: }
730: }
732: static void FormMomentumRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[],PetscScalar fz[])
733: {
734: PetscInt ngp;
735: PetscScalar gp_xi[GAUSS_POINTS][NSD];
736: PetscScalar gp_weight[GAUSS_POINTS];
737: PetscInt p,i;
738: PetscScalar Ni_p[NODES_PER_EL];
739: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
740: PetscScalar J_p,fac;
742: /* define quadrature rule */
743: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
745: /* evaluate integral */
746: for (p = 0; p < ngp; p++) {
747: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
748: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
749: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
750: fac = gp_weight[p]*J_p;
752: for (i = 0; i < NODES_PER_EL; i++) {
753: Fe[NSD*i] -= fac*Ni_p[i]*fx[p];
754: Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
755: Fe[NSD*i+2] -= fac*Ni_p[i]*fz[p];
756: }
757: }
758: }
760: static void FormContinuityRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar hc[])
761: {
762: PetscInt ngp;
763: PetscScalar gp_xi[GAUSS_POINTS][NSD];
764: PetscScalar gp_weight[GAUSS_POINTS];
765: PetscInt p,i;
766: PetscScalar Ni_p[NODES_PER_EL];
767: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
768: PetscScalar J_p,fac;
770: /* define quadrature rule */
771: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
773: /* evaluate integral */
774: for (p = 0; p < ngp; p++) {
775: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
776: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
777: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
778: fac = gp_weight[p]*J_p;
780: for (i = 0; i < NODES_PER_EL; i++) Fe[i] -= fac*Ni_p[i]*hc[p];
781: }
782: }
784: #define _ZERO_ROWCOL_i(A,i) { \
785: PetscInt KK; \
786: PetscScalar tmp = A[24*(i)+(i)]; \
787: for (KK=0;KK<24;KK++) A[24*(i)+KK]=0.0; \
788: for (KK=0;KK<24;KK++) A[24*KK+(i)]=0.0; \
789: A[24*(i)+(i)] = tmp;} \
791: #define _ZERO_ROW_i(A,i) { \
792: PetscInt KK; \
793: for (KK=0;KK<8;KK++) A[8*(i)+KK]=0.0;}
795: #define _ZERO_COL_i(A,i) { \
796: PetscInt KK; \
797: for (KK=0;KK<8;KK++) A[24*KK+(i)]=0.0;}
799: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,CellProperties cell_properties)
800: {
801: DM cda;
802: Vec coords;
803: DMDACoor3d ***_coords;
804: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
805: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
806: PetscInt sex,sey,sez,mx,my,mz;
807: PetscInt ei,ej,ek;
808: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
809: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
810: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
811: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
812: PetscScalar el_coords[NODES_PER_EL*NSD];
813: GaussPointCoefficients *props;
814: PetscScalar *prop_eta;
815: PetscInt n,M,N,P;
816: PetscErrorCode ierr;
819: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
820: /* setup for coords */
821: DMGetCoordinateDM(stokes_da,&cda);
822: DMGetCoordinatesLocal(stokes_da,&coords);
823: DMDAVecGetArray(cda,coords,&_coords);
825: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
826: for (ek = sez; ek < sez+mz; ek++) {
827: for (ej = sey; ej < sey+my; ej++) {
828: for (ei = sex; ei < sex+mx; ei++) {
829: /* get coords for the element */
830: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
831: /* get cell properties */
832: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
833: /* get coefficients for the element */
834: prop_eta = props->eta;
836: /* initialise element stiffness matrix */
837: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
838: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
839: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
840: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
842: /* form element stiffness matrix */
843: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
844: FormGradientOperatorQ13D(Ge,el_coords);
845: /*#if defined(ASSEMBLE_LOWER_TRIANGULAR)*/
846: FormDivergenceOperatorQ13D(De,el_coords);
847: /*#endif*/
848: FormStabilisationOperatorQ13D(Ce,el_coords,prop_eta);
850: /* insert element matrix into global matrix */
851: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
853: for (n=0; n<NODES_PER_EL; n++) {
854: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
855: _ZERO_ROWCOL_i(Ae,3*n);
856: _ZERO_ROW_i (Ge,3*n);
857: _ZERO_COL_i (De,3*n);
858: }
860: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
861: _ZERO_ROWCOL_i(Ae,3*n+1);
862: _ZERO_ROW_i (Ge,3*n+1);
863: _ZERO_COL_i (De,3*n+1);
864: }
866: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
867: _ZERO_ROWCOL_i(Ae,3*n+2);
868: _ZERO_ROW_i (Ge,3*n+2);
869: _ZERO_COL_i (De,3*n+2);
870: }
871: }
872: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
873: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
874: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
875: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
876: }
877: }
878: }
879: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
880: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
882: DMDAVecRestoreArray(cda,coords,&_coords);
884: return(0);
885: }
887: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,CellProperties cell_properties)
888: {
889: DM cda;
890: Vec coords;
891: DMDACoor3d ***_coords;
892: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
893: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
894: PetscInt sex,sey,sez,mx,my,mz;
895: PetscInt ei,ej,ek;
896: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
897: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
898: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
899: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
900: PetscScalar el_coords[NODES_PER_EL*NSD];
901: GaussPointCoefficients *props;
902: PetscScalar *prop_eta;
903: PetscInt n,M,N,P;
904: PetscErrorCode ierr;
907: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
908: /* setup for coords */
909: DMGetCoordinateDM(stokes_da,&cda);
910: DMGetCoordinatesLocal(stokes_da,&coords);
911: DMDAVecGetArray(cda,coords,&_coords);
913: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
914: for (ek = sez; ek < sez+mz; ek++) {
915: for (ej = sey; ej < sey+my; ej++) {
916: for (ei = sex; ei < sex+mx; ei++) {
917: /* get coords for the element */
918: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
919: /* get cell properties */
920: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
921: /* get coefficients for the element */
922: prop_eta = props->eta;
924: /* initialise element stiffness matrix */
925: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
926: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
927: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
928: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
930: /* form element stiffness matrix */
931: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
932: FormGradientOperatorQ13D(Ge,el_coords);
933: /* FormDivergenceOperatorQ13D(De,el_coords); */
934: FormScaledMassMatrixOperatorQ13D(Ce,el_coords,prop_eta);
936: /* insert element matrix into global matrix */
937: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
939: for (n=0; n<NODES_PER_EL; n++) {
940: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
941: _ZERO_ROWCOL_i(Ae,3*n);
942: _ZERO_ROW_i (Ge,3*n);
943: _ZERO_COL_i (De,3*n);
944: }
946: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
947: _ZERO_ROWCOL_i(Ae,3*n+1);
948: _ZERO_ROW_i (Ge,3*n+1);
949: _ZERO_COL_i (De,3*n+1);
950: }
952: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
953: _ZERO_ROWCOL_i(Ae,3*n+2);
954: _ZERO_ROW_i (Ge,3*n+2);
955: _ZERO_COL_i (De,3*n+2);
956: }
957: }
958: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
959: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
960: /*MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);*/
961: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
962: }
963: }
964: }
965: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
966: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
968: DMDAVecRestoreArray(cda,coords,&_coords);
969: return(0);
970: }
972: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,CellProperties cell_properties)
973: {
974: DM cda;
975: Vec coords;
976: DMDACoor3d ***_coords;
977: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
978: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
979: PetscInt sex,sey,sez,mx,my,mz;
980: PetscInt ei,ej,ek;
981: PetscScalar Fe[NODES_PER_EL*U_DOFS];
982: PetscScalar He[NODES_PER_EL*P_DOFS];
983: PetscScalar el_coords[NODES_PER_EL*NSD];
984: GaussPointCoefficients *props;
985: PetscScalar *prop_fx,*prop_fy,*prop_fz,*prop_hc;
986: Vec local_F;
987: StokesDOF ***ff;
988: PetscInt n,M,N,P;
989: PetscErrorCode ierr;
992: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
993: /* setup for coords */
994: DMGetCoordinateDM(stokes_da,&cda);
995: DMGetCoordinatesLocal(stokes_da,&coords);
996: DMDAVecGetArray(cda,coords,&_coords);
998: /* get acces to the vector */
999: DMGetLocalVector(stokes_da,&local_F);
1000: VecZeroEntries(local_F);
1001: DMDAVecGetArray(stokes_da,local_F,&ff);
1003: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
1004: for (ek = sez; ek < sez+mz; ek++) {
1005: for (ej = sey; ej < sey+my; ej++) {
1006: for (ei = sex; ei < sex+mx; ei++) {
1007: /* get coords for the element */
1008: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1009: /* get cell properties */
1010: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
1011: /* get coefficients for the element */
1012: prop_fx = props->fx;
1013: prop_fy = props->fy;
1014: prop_fz = props->fz;
1015: prop_hc = props->hc;
1017: /* initialise element stiffness matrix */
1018: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
1019: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
1021: /* form element stiffness matrix */
1022: FormMomentumRhsQ13D(Fe,el_coords,prop_fx,prop_fy,prop_fz);
1023: FormContinuityRhsQ13D(He,el_coords,prop_hc);
1025: /* insert element matrix into global matrix */
1026: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
1028: for (n=0; n<NODES_PER_EL; n++) {
1029: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) Fe[3*n] = 0.0;
1031: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) Fe[3*n+1] = 0.0;
1033: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) Fe[3*n+2] = 0.0;
1034: }
1036: DMDASetValuesLocalStencil3D_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
1037: }
1038: }
1039: }
1040: DMDAVecRestoreArray(stokes_da,local_F,&ff);
1041: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
1042: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
1043: DMRestoreLocalVector(stokes_da,&local_F);
1045: DMDAVecRestoreArray(cda,coords,&_coords);
1046: return(0);
1047: }
1049: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta,PetscReal *MX,PetscReal *MY,PetscReal *MZ)
1050: {
1051: *theta = 0.0;
1052: *MX = 2.0 * PETSC_PI;
1053: *MY = 2.0 * PETSC_PI;
1054: *MZ = 2.0 * PETSC_PI;
1055: }
1056: static void evaluate_MS_FrankKamentski(PetscReal pos[],PetscReal v[],PetscReal *p,PetscReal *eta,PetscReal Fm[],PetscReal *Fc)
1057: {
1058: PetscReal x,y,z;
1059: PetscReal theta,MX,MY,MZ;
1061: evaluate_MS_FrankKamentski_constants(&theta,&MX,&MY,&MZ);
1062: x = pos[0];
1063: y = pos[1];
1064: z = pos[2];
1065: if (v) {
1066: /*
1067: v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1068: v[1] = z*cos(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1069: v[2] = -(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z);
1070: */
1071: /*
1072: v[0] = PetscSinReal(PETSC_PI*x);
1073: v[1] = PetscSinReal(PETSC_PI*y);
1074: v[2] = PetscSinReal(PETSC_PI*z);
1075: */
1076: v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1077: v[1] = z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1078: v[2] = PetscPowRealInt(z,2)*(PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)/2 - PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)/2) - PETSC_PI*PetscPowRealInt(z,4)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y)/4;
1079: }
1080: if (p) *p = PetscPowRealInt(x,2) + PetscPowRealInt(y,2) + PetscPowRealInt(z,2);
1081: if (eta) {
1082: /**eta = PetscExpReal(-theta*(1.0 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));*/
1083: *eta = 1.0;
1084: }
1085: if (Fm) {
1086: /*
1087: Fm[0] = -2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(PETSC_PI*x) - 0.2*MZ*theta*(-1.5*PetscPowRealInt(x,2)*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PetscPowRealInt(z,2)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x))*PetscCosReal(MX*x)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) - 0.2*PETSC_PI*MX*theta*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscCosReal(MZ*z)*PetscExpReal(y)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) + 2.0*(3.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 3.0*PETSC_PI*PetscPowRealInt(x,2)*PetscCosReal(2.0*PETSC_PI*z))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 1.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) ;
1088: Fm[1] = -2*y - 0.2*MX*theta*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x))*PetscCosReal(MZ*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) - 0.2*MZ*theta*(-1.5*PetscPowRealInt(y,2)*PetscSinReal(2.0*PETSC_PI*z) + 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y))*PetscCosReal(MX*x)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) + 2.0*(-2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.5*PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 2*PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(-z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) - 6.0*PETSC_PI*PetscPowRealInt(y,2)*PetscCosReal(2.0*PETSC_PI*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));
1089: Fm[2] = -2*z + 8.0*PetscPowRealInt(PETSC_PI,2)*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(2.0*PETSC_PI*z) - 0.2*MX*theta*(-1.5*PetscPowRealInt(x,2)*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PetscPowRealInt(z,2)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x))*PetscCosReal(MZ*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) + 0.4*PETSC_PI*MZ*theta*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(MX*x)*PetscCosReal(2.0*PETSC_PI*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) + 2.0*(-3.0*x*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(-3.0*y*PetscSinReal(2.0*PETSC_PI*z) - 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.5*PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(-1.5*PetscPowRealInt(y,2)*PetscSinReal(2.0*PETSC_PI*z) + 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) ;
1090: */
1091: /*
1092: Fm[0]=-2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*x);
1093: Fm[1]=-2*y - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*y);
1094: Fm[2]=-2*z - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*z);
1095: */
1096: /*
1097: Fm[0] = -2*x + PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 6.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 6.0*PETSC_PI*PetscPowRealInt(x,2)*PetscCosReal(2.0*PETSC_PI*z) - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 2.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x) ;
1098: Fm[1] = -2*y - 6.0*PETSC_PI*PetscPowRealInt(y,2)*PetscCosReal(2.0*PETSC_PI*z) + 2.0*z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 6.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) - 4.0*PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);
1099: Fm[2] = -2*z - 6.0*x*PetscSinReal(2.0*PETSC_PI*z) - 6.0*y*PetscSinReal(2.0*PETSC_PI*z) - PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 8.0*PetscPowRealInt(PETSC_PI,2)*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z) + 3.0*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y) ;
1100: */
1101: Fm[0] = -2*x + 2*z*(PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x) - 1.0*PETSC_PI*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x)) + PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 6.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 2.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x);
1102: Fm[1] = -2*y + 2*z*(-PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)/2 + PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)/2 + PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)) + 2.0*z*PetscCosReal(2.0*PETSC_PI *x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 6.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 4.0*PETSC_PI*z*PetscCosReal(PETSC_PI *y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);
1103: Fm[2] = -2*z + PetscPowRealInt(z,2)*(-2.0*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 2.0*PetscPowRealInt(PETSC_PI,3)*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)) + PetscPowRealInt(z,2)*(PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)/2 - 3*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)/2 + PetscPowRealInt(PETSC_PI,3)*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)/2 - 3*PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)/2) + 1.0*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.25*PetscPowRealInt(PETSC_PI,3)*PetscPowRealInt(z,4)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) - 0.25*PETSC_PI*PetscPowRealInt(z,4)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) - 3.0*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) - 1.0*PETSC_PI*PetscCosReal(PETSC_PI *y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);
1104: }
1105: if (Fc) {
1106: /**Fc = -2.0*PETSC_PI*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(2.0*PETSC_PI*z) - z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y) ;*/
1107: /**Fc = PETSC_PI*PetscCosReal(PETSC_PI*x) + PETSC_PI*PetscCosReal(PETSC_PI*y) + PETSC_PI*PetscCosReal(PETSC_PI*z);*/
1108: /**Fc = -2.0*PETSC_PI*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(2.0*PETSC_PI*z) - z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);*/
1109: *Fc = 0.0;
1110: }
1111: }
1113: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx,PetscInt my,PetscInt mz,DM *_da,Vec *_X)
1114: {
1115: DM da,cda;
1116: Vec X;
1117: StokesDOF ***_stokes;
1118: Vec coords;
1119: DMDACoor3d ***_coords;
1120: PetscInt si,sj,sk,ei,ej,ek,i,j,k;
1124: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1125: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,4,1,NULL,NULL,NULL,&da);
1126: DMSetFromOptions(da);
1127: DMSetUp(da);
1128: DMDASetFieldName(da,0,"anlytic_Vx");
1129: DMDASetFieldName(da,1,"anlytic_Vy");
1130: DMDASetFieldName(da,2,"anlytic_Vz");
1131: DMDASetFieldName(da,3,"analytic_P");
1133: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
1135: DMGetCoordinatesLocal(da,&coords);
1136: DMGetCoordinateDM(da,&cda);
1137: DMDAVecGetArray(cda,coords,&_coords);
1139: DMCreateGlobalVector(da,&X);
1140: DMDAVecGetArray(da,X,&_stokes);
1142: DMDAGetCorners(da,&si,&sj,&sk,&ei,&ej,&ek);
1143: for (k = sk; k < sk+ek; k++) {
1144: for (j = sj; j < sj+ej; j++) {
1145: for (i = si; i < si+ei; i++) {
1146: PetscReal pos[NSD],pressure,vel[NSD];
1148: pos[0] = PetscRealPart(_coords[k][j][i].x);
1149: pos[1] = PetscRealPart(_coords[k][j][i].y);
1150: pos[2] = PetscRealPart(_coords[k][j][i].z);
1152: evaluate_MS_FrankKamentski(pos,vel,&pressure,NULL,NULL,NULL);
1154: _stokes[k][j][i].u_dof = vel[0];
1155: _stokes[k][j][i].v_dof = vel[1];
1156: _stokes[k][j][i].w_dof = vel[2];
1157: _stokes[k][j][i].p_dof = pressure;
1158: }
1159: }
1160: }
1161: DMDAVecRestoreArray(da,X,&_stokes);
1162: DMDAVecRestoreArray(cda,coords,&_coords);
1164: *_da = da;
1165: *_X = X;
1166: return(0);
1167: }
1169: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da,Vec X,Vec X_analytic)
1170: {
1171: DM cda;
1172: Vec coords,X_analytic_local,X_local;
1173: DMDACoor3d ***_coords;
1174: PetscInt sex,sey,sez,mx,my,mz;
1175: PetscInt ei,ej,ek;
1176: PetscScalar el_coords[NODES_PER_EL*NSD];
1177: StokesDOF ***stokes_analytic,***stokes;
1178: StokesDOF stokes_analytic_e[NODES_PER_EL],stokes_e[NODES_PER_EL];
1180: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1181: PetscScalar Ni_p[NODES_PER_EL];
1182: PetscInt ngp;
1183: PetscScalar gp_xi[GAUSS_POINTS][NSD];
1184: PetscScalar gp_weight[GAUSS_POINTS];
1185: PetscInt p,i;
1186: PetscScalar J_p,fac;
1187: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1188: PetscScalar tint_p_ms,tint_p,int_p_ms,int_p;
1189: PetscInt M;
1190: PetscReal xymin[NSD],xymax[NSD];
1194: /* define quadrature rule */
1195: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1197: /* setup for coords */
1198: DMGetCoordinateDM(stokes_da,&cda);
1199: DMGetCoordinatesLocal(stokes_da,&coords);
1200: DMDAVecGetArray(cda,coords,&_coords);
1202: /* setup for analytic */
1203: DMCreateLocalVector(stokes_da,&X_analytic_local);
1204: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1205: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1206: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1208: /* setup for solution */
1209: DMCreateLocalVector(stokes_da,&X_local);
1210: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1211: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1212: DMDAVecGetArray(stokes_da,X_local,&stokes);
1214: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1215: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1217: h = (xymax[0]-xymin[0])/((PetscReal)(M-1));
1219: tp_L2 = tu_L2 = tu_H1 = 0.0;
1220: tint_p_ms = tint_p = 0.0;
1222: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
1224: for (ek = sez; ek < sez+mz; ek++) {
1225: for (ej = sey; ej < sey+my; ej++) {
1226: for (ei = sex; ei < sex+mx; ei++) {
1227: /* get coords for the element */
1228: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1229: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1230: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1232: /* evaluate integral */
1233: for (p = 0; p < ngp; p++) {
1234: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1235: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1236: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1237: fac = gp_weight[p]*J_p;
1239: for (i = 0; i < NODES_PER_EL; i++) {
1240: tint_p_ms = tint_p_ms+fac*Ni_p[i]*stokes_analytic_e[i].p_dof;
1241: tint_p = tint_p +fac*Ni_p[i]*stokes_e[i].p_dof;
1242: }
1243: }
1244: }
1245: }
1246: }
1247: MPI_Allreduce(&tint_p_ms,&int_p_ms,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1248: MPI_Allreduce(&tint_p,&int_p,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1250: PetscPrintf(PETSC_COMM_WORLD,"\\int P dv %1.4e (h) %1.4e (ms)\n",PetscRealPart(int_p),PetscRealPart(int_p_ms));
1252: /* remove mine and add manufacture one */
1253: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1254: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1256: {
1257: PetscInt k,L,dof;
1258: PetscScalar *fields;
1259: DMDAGetInfo(stokes_da,0, 0,0,0, 0,0,0, &dof,0,0,0,0,0);
1261: VecGetLocalSize(X_local,&L);
1262: VecGetArray(X_local,&fields);
1263: for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1264: VecRestoreArray(X_local,&fields);
1266: VecGetLocalSize(X,&L);
1267: VecGetArray(X,&fields);
1268: for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1269: VecRestoreArray(X,&fields);
1270: }
1272: DMDAVecGetArray(stokes_da,X_local,&stokes);
1273: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1275: for (ek = sez; ek < sez+mz; ek++) {
1276: for (ej = sey; ej < sey+my; ej++) {
1277: for (ei = sex; ei < sex+mx; ei++) {
1278: /* get coords for the element */
1279: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1280: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1281: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1283: /* evaluate integral */
1284: p_e_L2 = 0.0;
1285: u_e_L2 = 0.0;
1286: u_e_H1 = 0.0;
1287: for (p = 0; p < ngp; p++) {
1288: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1289: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1290: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1291: fac = gp_weight[p]*J_p;
1293: for (i = 0; i < NODES_PER_EL; i++) {
1294: PetscScalar u_error,v_error,w_error;
1296: p_e_L2 = p_e_L2+fac*Ni_p[i]*(stokes_e[i].p_dof-stokes_analytic_e[i].p_dof)*(stokes_e[i].p_dof-stokes_analytic_e[i].p_dof);
1298: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1299: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1300: w_error = stokes_e[i].w_dof-stokes_analytic_e[i].w_dof;
1301: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error+w_error*w_error);
1303: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1304: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error
1305: +GNx_p[2][i]*u_error*GNx_p[2][i]*u_error
1306: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1307: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error
1308: +GNx_p[2][i]*v_error*GNx_p[2][i]*v_error
1309: +GNx_p[0][i]*w_error*GNx_p[0][i]*w_error /* dw/dx */
1310: +GNx_p[1][i]*w_error*GNx_p[1][i]*w_error
1311: +GNx_p[2][i]*w_error*GNx_p[2][i]*w_error);
1312: }
1313: }
1315: tp_L2 += p_e_L2;
1316: tu_L2 += u_e_L2;
1317: tu_H1 += u_e_H1;
1318: }
1319: }
1320: }
1321: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1322: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1323: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1324: p_L2 = PetscSqrtScalar(p_L2);
1325: u_L2 = PetscSqrtScalar(u_L2);
1326: u_H1 = PetscSqrtScalar(u_H1);
1328: PetscPrintf(PETSC_COMM_WORLD,"%1.4e %1.4e %1.4e %1.4e \n",PetscRealPart(h),PetscRealPart(p_L2),PetscRealPart(u_L2),PetscRealPart(u_H1));
1331: DMDAVecRestoreArray(cda,coords,&_coords);
1333: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1334: VecDestroy(&X_analytic_local);
1335: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1336: VecDestroy(&X_local);
1337: return(0);
1338: }
1340: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da,Vec FIELD,const char file_prefix[])
1341: {
1342: char vtk_filename[PETSC_MAX_PATH_LEN];
1343: PetscMPIInt rank;
1344: MPI_Comm comm;
1345: FILE *vtk_fp = NULL;
1346: PetscInt si,sj,sk,nx,ny,nz,i;
1347: PetscInt f,n_fields,N;
1348: DM cda;
1349: Vec coords;
1350: Vec l_FIELD;
1351: PetscScalar *_L_FIELD;
1352: PetscInt memory_offset;
1353: PetscScalar *buffer;
1358: /* create file name */
1359: PetscObjectGetComm((PetscObject)da,&comm);
1360: MPI_Comm_rank(comm,&rank);
1361: PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"subdomain-%s-p%1.4d.vts",file_prefix,rank);
1363: /* open file and write header */
1364: vtk_fp = fopen(vtk_filename,"w");
1365: if (!vtk_fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"Cannot open file = %s \n",vtk_filename);
1367: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");
1369: /* coords */
1370: DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1371: N = nx * ny * nz;
1373: #if defined(PETSC_WORDS_BIGENDIAN)
1374: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1375: #else
1376: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1377: #endif
1378: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",si,si+nx-1,sj,sj+ny-1,sk,sk+nz-1);
1379: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",si,si+nx-1,sj,sj+ny-1,sk,sk+nz-1);
1381: memory_offset = 0;
1383: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <CellData></CellData>\n");
1385: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <Points>\n");
1387: /* copy coordinates */
1388: DMGetCoordinateDM(da,&cda);
1389: DMGetCoordinatesLocal(da,&coords);
1390: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%d\" />\n",memory_offset);
1391: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N*3;
1393: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </Points>\n");
1395: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PointData Scalars=\" ");
1396: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_fields,0,0,0,0,0);
1397: for (f=0; f<n_fields; f++) {
1398: const char *field_name;
1399: DMDAGetFieldName(da,f,&field_name);
1400: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"%s ",field_name);
1401: }
1402: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\">\n");
1404: for (f=0; f<n_fields; f++) {
1405: const char *field_name;
1407: DMDAGetFieldName(da,f,&field_name);
1408: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <DataArray type=\"Float64\" Name=\"%s\" format=\"appended\" offset=\"%d\"/>\n", field_name,memory_offset);
1409: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N;
1410: }
1412: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PointData>\n");
1413: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </Piece>\n");
1414: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </StructuredGrid>\n");
1416: PetscMalloc1(N,&buffer);
1417: DMGetLocalVector(da,&l_FIELD);
1418: DMGlobalToLocalBegin(da, FIELD,INSERT_VALUES,l_FIELD);
1419: DMGlobalToLocalEnd(da,FIELD,INSERT_VALUES,l_FIELD);
1420: VecGetArray(l_FIELD,&_L_FIELD);
1422: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <AppendedData encoding=\"raw\">\n");
1423: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"_");
1425: /* write coordinates */
1426: {
1427: int length = sizeof(PetscScalar)*N*3;
1428: PetscScalar *allcoords;
1430: fwrite(&length,sizeof(int),1,vtk_fp);
1431: VecGetArray(coords,&allcoords);
1432: fwrite(allcoords,sizeof(PetscScalar),3*N,vtk_fp);
1433: VecRestoreArray(coords,&allcoords);
1434: }
1435: /* write fields */
1436: for (f=0; f<n_fields; f++) {
1437: int length = sizeof(PetscScalar)*N;
1438: fwrite(&length,sizeof(int),1,vtk_fp);
1439: /* load */
1440: for (i=0; i<N; i++) buffer[i] = _L_FIELD[n_fields*i + f];
1442: /* write */
1443: fwrite(buffer,sizeof(PetscScalar),N,vtk_fp);
1444: }
1445: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\n </AppendedData>\n");
1447: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");
1449: PetscFree(buffer);
1450: VecRestoreArray(l_FIELD,&_L_FIELD);
1451: DMRestoreLocalVector(da,&l_FIELD);
1453: if (vtk_fp) {
1454: fclose(vtk_fp);
1455: vtk_fp = NULL;
1456: }
1458: return(0);
1459: }
1461: PetscErrorCode DAViewVTK_write_PieceExtend(FILE *vtk_fp,PetscInt indent_level,DM da,const char local_file_prefix[])
1462: {
1463: PetscMPIInt size,rank;
1464: MPI_Comm comm;
1465: const PetscInt *lx,*ly,*lz;
1466: PetscInt M,N,P,pM,pN,pP,sum,*olx,*oly,*olz;
1467: PetscInt *osx,*osy,*osz,*oex,*oey,*oez;
1468: PetscInt i,j,k,II,stencil;
1472: /* create file name */
1473: PetscObjectGetComm((PetscObject)da,&comm);
1474: MPI_Comm_size(comm,&size);
1475: MPI_Comm_rank(comm,&rank);
1477: DMDAGetInfo(da,0,&M,&N,&P,&pM,&pN,&pP,0,&stencil,0,0,0,0);
1478: DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
1480: /* generate start,end list */
1481: PetscMalloc1(pM+1,&olx);
1482: PetscMalloc1(pN+1,&oly);
1483: PetscMalloc1(pP+1,&olz);
1484: sum = 0;
1485: for (i=0; i<pM; i++) {
1486: olx[i] = sum;
1487: sum = sum + lx[i];
1488: }
1489: olx[pM] = sum;
1490: sum = 0;
1491: for (i=0; i<pN; i++) {
1492: oly[i] = sum;
1493: sum = sum + ly[i];
1494: }
1495: oly[pN] = sum;
1496: sum = 0;
1497: for (i=0; i<pP; i++) {
1498: olz[i] = sum;
1499: sum = sum + lz[i];
1500: }
1501: olz[pP] = sum;
1503: PetscMalloc1(pM,&osx);
1504: PetscMalloc1(pN,&osy);
1505: PetscMalloc1(pP,&osz);
1506: PetscMalloc1(pM,&oex);
1507: PetscMalloc1(pN,&oey);
1508: PetscMalloc1(pP,&oez);
1509: for (i=0; i<pM; i++) {
1510: osx[i] = olx[i] - stencil;
1511: oex[i] = olx[i] + lx[i] + stencil;
1512: if (osx[i]<0) osx[i]=0;
1513: if (oex[i]>M) oex[i]=M;
1514: }
1516: for (i=0; i<pN; i++) {
1517: osy[i] = oly[i] - stencil;
1518: oey[i] = oly[i] + ly[i] + stencil;
1519: if (osy[i]<0)osy[i]=0;
1520: if (oey[i]>M)oey[i]=N;
1521: }
1522: for (i=0; i<pP; i++) {
1523: osz[i] = olz[i] - stencil;
1524: oez[i] = olz[i] + lz[i] + stencil;
1525: if (osz[i]<0) osz[i]=0;
1526: if (oez[i]>P) oez[i]=P;
1527: }
1529: for (k=0; k<pP; k++) {
1530: for (j=0; j<pN; j++) {
1531: for (i=0; i<pM; i++) {
1532: char name[PETSC_MAX_PATH_LEN];
1533: PetscInt procid = i + j*pM + k*pM*pN; /* convert proc(i,j,k) to pid */
1534: PetscSNPrintf(name,sizeof(name),"subdomain-%s-p%1.4d.vts",local_file_prefix,procid);
1535: for (II=0; II<indent_level; II++) PetscFPrintf(PETSC_COMM_SELF,vtk_fp," ");
1537: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<Piece Extent=\"%d %d %d %d %d %d\" Source=\"%s\"/>\n",
1538: osx[i],oex[i]-1,
1539: osy[j],oey[j]-1,
1540: osz[k],oez[k]-1,name);
1541: }
1542: }
1543: }
1544: PetscFree(olx);
1545: PetscFree(oly);
1546: PetscFree(olz);
1547: PetscFree(osx);
1548: PetscFree(osy);
1549: PetscFree(osz);
1550: PetscFree(oex);
1551: PetscFree(oey);
1552: PetscFree(oez);
1553: return(0);
1554: }
1556: PetscErrorCode DAView_3DVTK_PStructuredGrid(DM da,const char file_prefix[],const char local_file_prefix[])
1557: {
1558: MPI_Comm comm;
1559: PetscMPIInt size,rank;
1560: char vtk_filename[PETSC_MAX_PATH_LEN];
1561: FILE *vtk_fp = NULL;
1562: PetscInt M,N,P,si,sj,sk,nx,ny,nz;
1563: PetscInt i,dofs;
1567: /* only master generates this file */
1568: PetscObjectGetComm((PetscObject)da,&comm);
1569: MPI_Comm_size(comm,&size);
1570: MPI_Comm_rank(comm,&rank);
1572: if (rank != 0) return(0);
1574: /* create file name */
1575: PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"%s.pvts",file_prefix);
1576: vtk_fp = fopen(vtk_filename,"w");
1577: if (!vtk_fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"Cannot open file = %s \n",vtk_filename);
1579: /* (VTK) generate pvts header */
1580: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");
1582: #if defined(PETSC_WORDS_BIGENDIAN)
1583: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1584: #else
1585: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1586: #endif
1588: /* define size of the nodal mesh based on the cell DM */
1589: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,&dofs,0,0,0,0,0);
1590: DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1591: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PStructuredGrid GhostLevel=\"1\" WholeExtent=\"%d %d %d %d %d %d\">\n",0,M-1,0,N-1,0,P-1); /* note overlap = 1 for Q1 */
1593: /* DUMP THE CELL REFERENCES */
1594: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PCellData>\n");
1595: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PCellData>\n");
1597: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PPoints>\n");
1598: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PDataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"%d\"/>\n",NSD);
1599: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PPoints>\n");
1601: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PPointData>\n");
1602: for (i=0; i<dofs; i++) {
1603: const char *fieldname;
1604: DMDAGetFieldName(da,i,&fieldname);
1605: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n",fieldname);
1606: }
1607: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PPointData>\n");
1609: /* write out the parallel information */
1610: DAViewVTK_write_PieceExtend(vtk_fp,2,da,local_file_prefix);
1612: /* close the file */
1613: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PStructuredGrid>\n");
1614: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");
1616: if (vtk_fp) {
1617: fclose(vtk_fp);
1618: vtk_fp = NULL;
1619: }
1620: return(0);
1621: }
1623: PetscErrorCode DAView3DPVTS(DM da, Vec x,const char NAME[])
1624: {
1625: char vts_filename[PETSC_MAX_PATH_LEN];
1626: char pvts_filename[PETSC_MAX_PATH_LEN];
1630: PetscSNPrintf(vts_filename,sizeof(vts_filename),"%s-mesh",NAME);
1631: DAView_3DVTK_StructuredGrid_appended(da,x,vts_filename);
1633: PetscSNPrintf(pvts_filename,sizeof(pvts_filename),"%s-mesh",NAME);
1634: DAView_3DVTK_PStructuredGrid(da,pvts_filename,vts_filename);
1635: return(0);
1636: }
1638: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
1639: {
1641: PetscReal norms[4];
1642: Vec Br,v,w;
1643: Mat A;
1646: KSPGetOperators(ksp,&A,NULL);
1647: MatCreateVecs(A,&w,&v);
1649: KSPBuildResidual(ksp,v,w,&Br);
1651: VecStrideNorm(Br,0,NORM_2,&norms[0]);
1652: VecStrideNorm(Br,1,NORM_2,&norms[1]);
1653: VecStrideNorm(Br,2,NORM_2,&norms[2]);
1654: VecStrideNorm(Br,3,NORM_2,&norms[3]);
1656: VecDestroy(&v);
1657: VecDestroy(&w);
1659: PetscPrintf(PETSC_COMM_WORLD,"%3D KSP Component U,V,W,P residual norm [ %1.12e, %1.12e, %1.12e, %1.12e ]\n",n,norms[0],norms[1],norms[2],norms[3]);
1660: return(0);
1661: }
1663: static PetscErrorCode PCMGSetupViaCoarsen(PC pc,DM da_fine)
1664: {
1665: PetscInt nlevels,k;
1666: PETSC_UNUSED PetscInt finest;
1667: DM *da_list,*daclist;
1668: Mat R;
1669: PetscErrorCode ierr;
1672: nlevels = 1;
1673: PetscOptionsGetInt(NULL,NULL,"-levels",&nlevels,0);
1675: PetscMalloc1(nlevels,&da_list);
1676: for (k=0; k<nlevels; k++) da_list[k] = NULL;
1677: PetscMalloc1(nlevels,&daclist);
1678: for (k=0; k<nlevels; k++) daclist[k] = NULL;
1680: /* finest grid is nlevels - 1 */
1681: finest = nlevels - 1;
1682: daclist[0] = da_fine;
1683: PetscObjectReference((PetscObject)da_fine);
1684: DMCoarsenHierarchy(da_fine,nlevels-1,&daclist[1]);
1685: for (k=0; k<nlevels; k++) {
1686: da_list[k] = daclist[nlevels-1-k];
1687: DMDASetUniformCoordinates(da_list[k],0.0,1.0,0.0,1.0,0.0,1.0);
1688: }
1690: PCMGSetLevels(pc,nlevels,NULL);
1691: PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
1692: PCMGSetGalerkin(pc,PC_MG_GALERKIN_PMAT);
1694: for (k=1; k<nlevels; k++) {
1695: DMCreateInterpolation(da_list[k-1],da_list[k],&R,NULL);
1696: PCMGSetInterpolation(pc,k,R);
1697: MatDestroy(&R);
1698: }
1700: /* tidy up */
1701: for (k=0; k<nlevels; k++) {
1702: DMDestroy(&da_list[k]);
1703: }
1704: PetscFree(da_list);
1705: PetscFree(daclist);
1706: return(0);
1707: }
1709: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx,PetscInt my,PetscInt mz)
1710: {
1711: DM da_Stokes;
1712: PetscInt u_dof,p_dof,dof,stencil_width;
1713: Mat A,B;
1714: PetscInt mxl,myl,mzl;
1715: DM vel_cda;
1716: Vec vel_coords;
1717: PetscInt p;
1718: Vec f,X;
1719: DMDACoor3d ***_vel_coords;
1720: PetscInt its;
1721: KSP ksp_S;
1722: PetscInt model_definition = 0;
1723: PetscInt ei,ej,ek,sex,sey,sez,Imx,Jmy,Kmz;
1724: CellProperties cell_properties;
1725: PetscBool write_output = PETSC_FALSE;
1729: /* Generate the da for velocity and pressure */
1730: /* Num nodes in each direction is mx+1, my+1, mz+1 */
1731: u_dof = U_DOFS; /* Vx, Vy - velocities */
1732: p_dof = P_DOFS; /* p - pressure */
1733: dof = u_dof+p_dof;
1734: stencil_width = 1;
1735: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1736: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,NULL,&da_Stokes);
1737: DMSetFromOptions(da_Stokes);
1738: DMSetUp(da_Stokes);
1739: DMDASetFieldName(da_Stokes,0,"Vx");
1740: DMDASetFieldName(da_Stokes,1,"Vy");
1741: DMDASetFieldName(da_Stokes,2,"Vz");
1742: DMDASetFieldName(da_Stokes,3,"P");
1744: /* unit box [0,1] x [0,1] x [0,1] */
1745: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.0,1.0);
1747: /* local number of elements */
1748: DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,&mzl);
1750: /* create quadrature point info for PDE coefficients */
1751: CellPropertiesCreate(da_Stokes,&cell_properties);
1753: /* interpolate the coordinates to quadrature points */
1754: DMGetCoordinateDM(da_Stokes,&vel_cda);
1755: DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1756: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1757: DMDAGetElementCorners(da_Stokes,&sex,&sey,&sez,&Imx,&Jmy,&Kmz);
1758: for (ek = sez; ek < sez+Kmz; ek++) {
1759: for (ej = sey; ej < sey+Jmy; ej++) {
1760: for (ei = sex; ei < sex+Imx; ei++) {
1761: /* get coords for the element */
1762: PetscInt ngp;
1763: PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
1764: PetscScalar el_coords[NSD*NODES_PER_EL];
1765: GaussPointCoefficients *cell;
1767: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1768: GetElementCoords3D(_vel_coords,ei,ej,ek,el_coords);
1769: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1771: for (p = 0; p < GAUSS_POINTS; p++) {
1772: PetscScalar xi_p[NSD],Ni_p[NODES_PER_EL];
1773: PetscScalar gp_x,gp_y,gp_z;
1774: PetscInt n;
1776: xi_p[0] = gp_xi[p][0];
1777: xi_p[1] = gp_xi[p][1];
1778: xi_p[2] = gp_xi[p][2];
1779: ShapeFunctionQ13D_Evaluate(xi_p,Ni_p);
1781: gp_x = gp_y = gp_z = 0.0;
1782: for (n = 0; n < NODES_PER_EL; n++) {
1783: gp_x = gp_x+Ni_p[n]*el_coords[NSD*n];
1784: gp_y = gp_y+Ni_p[n]*el_coords[NSD*n+1];
1785: gp_z = gp_z+Ni_p[n]*el_coords[NSD*n+2];
1786: }
1787: cell->gp_coords[NSD*p] = gp_x;
1788: cell->gp_coords[NSD*p+1] = gp_y;
1789: cell->gp_coords[NSD*p+2] = gp_z;
1790: }
1791: }
1792: }
1793: }
1795: PetscOptionsGetInt(NULL,NULL,"-model",&model_definition,NULL);
1797: switch (model_definition) {
1798: case 0: /* isoviscous */
1799: for (ek = sez; ek < sez+Kmz; ek++) {
1800: for (ej = sey; ej < sey+Jmy; ej++) {
1801: for (ei = sex; ei < sex+Imx; ei++) {
1802: GaussPointCoefficients *cell;
1804: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1805: for (p = 0; p < GAUSS_POINTS; p++) {
1806: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1807: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1808: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1810: cell->eta[p] = 1.0;
1812: cell->fx[p] = 0.0*coord_x;
1813: cell->fy[p] = -PetscSinReal(2.2*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1814: cell->fz[p] = 0.0*coord_z;
1815: cell->hc[p] = 0.0;
1816: }
1817: }
1818: }
1819: }
1820: break;
1822: case 1: /* manufactured */
1823: for (ek = sez; ek < sez+Kmz; ek++) {
1824: for (ej = sey; ej < sey+Jmy; ej++) {
1825: for (ei = sex; ei < sex+Imx; ei++) {
1826: PetscReal eta,Fm[NSD],Fc,pos[NSD];
1827: GaussPointCoefficients *cell;
1829: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1830: for (p = 0; p < GAUSS_POINTS; p++) {
1831: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1832: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1833: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1835: pos[0] = coord_x;
1836: pos[1] = coord_y;
1837: pos[2] = coord_z;
1839: evaluate_MS_FrankKamentski(pos,NULL,NULL,&eta,Fm,&Fc);
1840: cell->eta[p] = eta;
1842: cell->fx[p] = Fm[0];
1843: cell->fy[p] = Fm[1];
1844: cell->fz[p] = Fm[2];
1845: cell->hc[p] = Fc;
1846: }
1847: }
1848: }
1849: }
1850: break;
1852: case 2: /* solcx */
1853: for (ek = sez; ek < sez+Kmz; ek++) {
1854: for (ej = sey; ej < sey+Jmy; ej++) {
1855: for (ei = sex; ei < sex+Imx; ei++) {
1856: GaussPointCoefficients *cell;
1858: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1859: for (p = 0; p < GAUSS_POINTS; p++) {
1860: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1861: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1862: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1864: cell->eta[p] = 1.0;
1866: cell->fx[p] = 0.0;
1867: cell->fy[p] = -PetscSinReal(3.0*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1868: cell->fz[p] = 0.0*coord_z;
1869: cell->hc[p] = 0.0;
1870: }
1871: }
1872: }
1873: }
1874: break;
1876: case 3: /* sinker */
1877: for (ek = sez; ek < sez+Kmz; ek++) {
1878: for (ej = sey; ej < sey+Jmy; ej++) {
1879: for (ei = sex; ei < sex+Imx; ei++) {
1880: GaussPointCoefficients *cell;
1882: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1883: for (p = 0; p < GAUSS_POINTS; p++) {
1884: PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p]);
1885: PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);
1886: PetscReal zp = PetscRealPart(cell->gp_coords[NSD*p+2]);
1888: cell->eta[p] = 1.0e-2;
1889: cell->fx[p] = 0.0;
1890: cell->fy[p] = 0.0;
1891: cell->fz[p] = 0.0;
1892: cell->hc[p] = 0.0;
1894: if ((PetscAbsReal(xp-0.5) < 0.2) &&
1895: (PetscAbsReal(yp-0.5) < 0.2) &&
1896: (PetscAbsReal(zp-0.5) < 0.2)) {
1897: cell->eta[p] = 1.0;
1898: cell->fz[p] = 1.0;
1899: }
1901: }
1902: }
1903: }
1904: }
1905: break;
1907: default:
1908: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"No default model is supported. Choose either -model {0,1,2,3}");
1909: break;
1910: }
1912: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1914: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1915: DMSetMatType(da_Stokes,MATAIJ);
1916: DMCreateMatrix(da_Stokes,&A);
1917: DMCreateMatrix(da_Stokes,&B);
1918: DMCreateGlobalVector(da_Stokes,&X);
1919: DMCreateGlobalVector(da_Stokes,&f);
1921: /* assemble A11 */
1922: MatZeroEntries(A);
1923: MatZeroEntries(B);
1924: VecZeroEntries(f);
1926: AssembleA_Stokes(A,da_Stokes,cell_properties);
1927: AssembleA_PCStokes(B,da_Stokes,cell_properties);
1928: /* build force vector */
1929: AssembleF_Stokes(f,da_Stokes,cell_properties);
1931: /* SOLVE */
1932: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1933: KSPSetOptionsPrefix(ksp_S,"stokes_"); /* stokes */
1934: KSPSetOperators(ksp_S,A,B);
1935: KSPSetFromOptions(ksp_S);
1937: {
1938: PC pc;
1939: const PetscInt ufields[] = {0,1,2},pfields[] = {3};
1940: KSPGetPC(ksp_S,&pc);
1941: PCFieldSplitSetBlockSize(pc,4);
1942: PCFieldSplitSetFields(pc,"u",3,ufields,ufields);
1943: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1944: }
1946: {
1947: PC pc;
1948: PetscBool same = PETSC_FALSE;
1949: KSPGetPC(ksp_S,&pc);
1950: PetscObjectTypeCompare((PetscObject)pc,PCMG,&same);
1951: if (same) {
1952: PCMGSetupViaCoarsen(pc,da_Stokes);
1953: }
1954: }
1956: {
1957: PetscBool stokes_monitor = PETSC_FALSE;
1958: PetscOptionsGetBool(NULL,NULL,"-stokes_ksp_monitor_blocks",&stokes_monitor,0);
1959: if (stokes_monitor) {
1960: KSPMonitorSet(ksp_S,KSPMonitorStokesBlocks,NULL,NULL);
1961: }
1962: }
1963: KSPSolve(ksp_S,f,X);
1965: PetscOptionsGetBool(NULL,NULL,"-write_pvts",&write_output,NULL);
1966: if (write_output) {DAView3DPVTS(da_Stokes,X,"up");}
1967: {
1968: PetscBool flg = PETSC_FALSE;
1969: char filename[PETSC_MAX_PATH_LEN];
1970: PetscOptionsGetString(NULL,NULL,"-write_binary",filename,sizeof(filename),&flg);
1971: if (flg) {
1972: PetscViewer viewer;
1973: /* PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer); */
1974: PetscViewerVTKOpen(PETSC_COMM_WORLD,"ex42.vts",FILE_MODE_WRITE,&viewer);
1975: VecView(X,viewer);
1976: PetscViewerDestroy(&viewer);
1977: }
1978: }
1979: KSPGetIterationNumber(ksp_S,&its);
1981: /* verify */
1982: if (model_definition == 1) {
1983: DM da_Stokes_analytic;
1984: Vec X_analytic;
1986: DMDACreateManufacturedSolution(mx,my,mz,&da_Stokes_analytic,&X_analytic);
1987: if (write_output) {
1988: DAView3DPVTS(da_Stokes_analytic,X_analytic,"ms");
1989: }
1990: DMDAIntegrateErrors3D(da_Stokes_analytic,X,X_analytic);
1991: if (write_output) {
1992: DAView3DPVTS(da_Stokes,X,"up2");
1993: }
1994: DMDestroy(&da_Stokes_analytic);
1995: VecDestroy(&X_analytic);
1996: }
1998: KSPDestroy(&ksp_S);
1999: VecDestroy(&X);
2000: VecDestroy(&f);
2001: MatDestroy(&A);
2002: MatDestroy(&B);
2004: CellPropertiesDestroy(&cell_properties);
2005: DMDestroy(&da_Stokes);
2006: return(0);
2007: }
2009: int main(int argc,char **args)
2010: {
2012: PetscInt mx,my,mz;
2014: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
2015: mx = my = mz = 10;
2016: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
2017: my = mx; mz = mx;
2018: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
2019: PetscOptionsGetInt(NULL,NULL,"-mz",&mz,NULL);
2020: solve_stokes_3d_coupled(mx,my,mz);
2021: PetscFinalize();
2022: return ierr;
2023: }