Actual source code: ex55.c
petsc-3.5.4 2015-05-23
1: static char help[] = "Allen-Cahn-2d problem for constant mobility and triangular elements.\n\
2: Runtime options include:\n\
3: -xmin <xmin>\n\
4: -xmax <xmax>\n\
5: -ymin <ymin>\n\
6: -T <T>, where <T> is the end time for the time domain simulation\n\
7: -dt <dt>,where <dt> is the step size for the numerical integration\n\
8: -gamma <gamma>\n\
9: -theta_c <theta_c>\n\n";
11: /*
12: Solves the linear system using a Schur complement solver based on PCLSC, solve the A00 block with hypre BoomerAMG
14: ./ex55 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_detect_saddle_point -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self -fieldsplit_1_ksp_type fgmres -fieldsplit_1_pc_type lsc -snes_vi_monitor -ksp_monitor_true_residual -fieldsplit_ksp_monitor -fieldsplit_0_pc_type hypre -da_grid_x 65 -da_grid_y 65 -snes_atol 1.e-11 -ksp_rtol 1.e-8
16: Solves the linear systems with multigrid on the entire system using a Schur complement based smoother (which is handled by PCFIELDSPLIT)
18: ./ex55 -ksp_type fgmres -pc_type mg -mg_levels_ksp_type fgmres -mg_levels_pc_type fieldsplit -mg_levels_pc_fieldsplit_detect_saddle_point -mg_levels_pc_fieldsplit_type schur -mg_levels_pc_fieldsplit_factorization_type full -mg_levels_pc_fieldsplit_schur_precondition user -mg_levels_fieldsplit_1_ksp_type gmres -mg_levels_fieldsplit_1_pc_type none -mg_levels_fieldsplit_0_ksp_type preonly -mg_levels_fieldsplit_0_pc_type sor -mg_levels_fieldsplit_0_pc_sor_forward -snes_vi_monitor -ksp_monitor_true_residual -pc_mg_levels 5 -pc_mg_galerkin -mg_levels_ksp_monitor -mg_levels_fieldsplit_ksp_monitor -mg_levels_ksp_max_it 2 -mg_levels_fieldsplit_ksp_max_it 5 -snes_atol 1.e-11 -mg_coarse_ksp_type preonly -mg_coarse_pc_type svd -da_grid_x 65 -da_grid_y 65 -ksp_rtol 1.e-8
20: */
22: #include <petscsnes.h>
23: #include <petscdmda.h>
25: typedef struct {
26: PetscReal dt,T; /* Time step and end time */
27: DM da;
28: Mat M; /* Jacobian matrix */
29: Mat M_0;
30: Vec q,u1,u2,u3,work1,work2,work3,work4;
31: PetscScalar epsilon; /* physics parameters */
32: PetscReal xmin,xmax,ymin,ymax;
33: } AppCtx;
35: PetscErrorCode GetParams(AppCtx*);
36: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
37: PetscErrorCode SetUpMatrices(AppCtx*);
38: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
39: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
40: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
41: PetscErrorCode Update_q(Vec,Vec,Vec,Vec,Mat,AppCtx*);
42: PetscErrorCode Update_u(Vec,Vec,Vec,Vec);
46: int main(int argc, char **argv)
47: {
49: Vec x,r; /* Solution and residual vectors */
50: SNES snes; /* Nonlinear solver context */
51: AppCtx user; /* Application context */
52: Vec xl,xu; /* Upper and lower bounds on variables */
53: Mat J;
54: PetscScalar t=0.0;
55: PetscViewer view_out, view_q, view1;
57: PetscInitialize(&argc,&argv, (char*)0, help);
59: /* Get physics and time parameters */
60: GetParams(&user);
61: /* Create a 2D DA with dof = 2 */
62: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,NULL,NULL,&user.da);
63: /* Set Element type (triangular) */
64: DMDASetElementType(user.da,DMDA_ELEMENT_P1);
66: /* Set x and y coordinates */
67: DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0.0,1.0);
69: /* Get global vector x from DM and duplicate vectors r,xl,xu */
70: DMCreateGlobalVector(user.da,&x);
71: VecDuplicate(x,&r);
72: VecDuplicate(x,&xl);
73: VecDuplicate(x,&xu);
74: VecDuplicate(x,&user.q);
76: /* Get Jacobian matrix structure from the da */
77: DMSetMatType(user.da,MATAIJ);
78: DMCreateMatrix(user.da,&user.M);
79: /* Form the jacobian matrix and M_0 */
80: SetUpMatrices(&user);
81: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
83: /* Create nonlinear solver context */
84: SNESCreate(PETSC_COMM_WORLD,&snes);
85: SNESSetDM(snes,user.da);
87: /* Set Function evaluation and jacobian evaluation routines */
88: SNESSetFunction(snes,r,FormFunction,(void*)&user);
89: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
91: /* Set the boundary conditions */
92: SetVariableBounds(user.da,xl,xu);
93: SNESVISetVariableBounds(snes,xl,xu);
94: SNESSetFromOptions(snes);
96: SetInitialGuess(x,&user);
97: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
98: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
99: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file1",FILE_MODE_WRITE,&view1);
100: /* Begin time loop */
101: while (t < user.T) {
102: VecView(user.u1,view1);
103: VecView(user.u2,view1);
104: VecView(user.u3,view1);
106: Update_q(user.q,user.u1,user.u2,user.u3,user.M_0,&user);
107: VecView(user.q,view_q);
108: SNESSolve(snes,NULL,x);
109: PetscInt its;
110: SNESGetIterationNumber(snes,&its);
111: PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4f in %d iterations\n",t,its);
112: Update_u(user.u1,user.u2,user.u3,x);
113: t = t + user.dt;
114: VecView(user.u1,view_out);
115: VecView(user.u2,view_out);
116: VecView(user.u3,view_out);
117: }
121: PetscViewerDestroy(&view_out);
122: PetscViewerDestroy(&view_q);
123: PetscViewerDestroy(&view1);
125: VecDestroy(&x);
126: VecDestroy(&r);
127: VecDestroy(&xl);
128: VecDestroy(&xu);
129: VecDestroy(&user.q);
130: VecDestroy(&user.u1);
131: VecDestroy(&user.u2);
132: VecDestroy(&user.u3);
133: VecDestroy(&user.work1);
134: VecDestroy(&user.work2);
135: VecDestroy(&user.work3);
136: VecDestroy(&user.work4);
137: MatDestroy(&user.M);
138: MatDestroy(&user.M_0);
139: MatDestroy(&J);
140: DMDestroy(&user.da);
141: SNESDestroy(&snes);
142: PetscFinalize();
143: return 0;
144: }
148: PetscErrorCode Update_u(Vec u1,Vec u2,Vec u3,Vec X)
149: {
151: PetscInt i,n;
152: PetscScalar *u1_arr,*u2_arr,*u3_arr,*x_arr;
155: VecGetLocalSize(u1,&n);
156: VecGetArray(u1,&u1_arr);
157: VecGetArray(u2,&u2_arr);
158: VecGetArray(u3,&u3_arr);
159: VecGetArray(X,&x_arr);
160: for (i=0; i<n; i++) {
161: u1_arr[i] = x_arr[4*i];
162: u2_arr[i] = x_arr[4*i+1];
163: u3_arr[i] = x_arr[4*i+2];
164: }
165: VecRestoreArray(u1,&u1_arr);
166: VecRestoreArray(u2,&u2_arr);
167: VecRestoreArray(u3,&u3_arr);
168: VecRestoreArray(X,&x_arr);
169: return(0);
170: }
174: PetscErrorCode Update_q(Vec q,Vec u1,Vec u2,Vec u3,Mat M_0,AppCtx *user)
175: {
177: PetscScalar *q_arr,*w_arr;
178: PetscInt i,n;
179: /*PetscViewer view_q;*/
182: VecSet(user->work1,user->dt/3);
183: /* VecView(user->work1,PETSC_VIEWER_STDOUT_WORLD);*/
184: MatMult(M_0,user->work1,user->work2);
185: /* VecView(user->work2,PETSC_VIEWER_STDOUT_WORLD);*/
187: MatMult(M_0,u1,user->work1);
188: MatMult(M_0,u1,user->work4);
189: /* VecView(u1,PETSC_VIEWER_STDOUT_WORLD);*/
190: /* VecView(user->work4,PETSC_VIEWER_STDOUT_WORLD);*/
191: VecScale(user->work1,-1.0-(user->dt));
192: VecAXPY(user->work1,1.0,user->work2);
194: VecGetLocalSize(u1,&n);
195: VecGetArray(q,&q_arr);
196: VecGetArray(user->work1,&w_arr);
197: for (i=0; i<n; i++) q_arr[4*i] = w_arr[i];
198: VecRestoreArray(q,&q_arr);
199: VecRestoreArray(user->work1,&w_arr);
201: MatMult(M_0,u2,user->work1);
202: VecScale(user->work1,-1.0-(user->dt));
203: VecAXPY(user->work1,1.0,user->work2);
205: VecGetArray(q,&q_arr);
206: VecGetArray(user->work1,&w_arr);
207: for (i=0; i<n; i++) q_arr[4*i+1] = w_arr[i];
208: VecRestoreArray(q,&q_arr);
209: VecRestoreArray(user->work1,&w_arr);
211: MatMult(M_0,u3,user->work1);
212: VecScale(user->work1,-1.0-(user->dt));
213: VecAXPY(user->work1,1.0,user->work2);
215: VecGetArray(q,&q_arr);
216: VecGetArray(user->work1,&w_arr);
217: for (i=0; i<n; i++) {
218: q_arr[4*i+2] = w_arr[i];
219: q_arr[4*i+3] = 1.0;
220: }
221: VecRestoreArray(q,&q_arr);
222: VecRestoreArray(user->work1,&w_arr);
223: return(0);
224: }
228: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
229: {
230: PetscErrorCode ierr;
231: PetscInt nele,nen,n,i;
232: const PetscInt *ele;
233: Vec coords, rand1, rand2;
234: const PetscScalar *_coords;
235: PetscScalar x[3],y[3];
236: PetscInt idx[3];
237: PetscScalar *xx,*w1,*w2,*u1,*u2,*u3;
238: PetscViewer view_out;
241: /* Get ghosted coordinates */
242: DMGetCoordinatesLocal(user->da,&coords);
243: VecDuplicate(user->u1,&rand1);
244: VecDuplicate(user->u1,&rand2);
245: VecSetRandom(rand1,NULL);
246: VecSetRandom(rand2,NULL);
248: VecGetLocalSize(X,&n);
249: VecGetArrayRead(coords,&_coords);
250: VecGetArray(X,&xx);
251: VecGetArray(user->work1,&w1);
252: VecGetArray(user->work2,&w2);
253: VecGetArray(user->u1,&u1);
254: VecGetArray(user->u2,&u2);
255: VecGetArray(user->u3,&u3);
257: /* Get local element info */
258: DMDAGetElements(user->da,&nele,&nen,&ele);
259: for (i=0; i < nele; i++) {
260: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
261: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
262: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
263: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
265: PetscScalar vals1[3],vals2[3],valsrand[3];
266: PetscInt r;
267: for (r=0; r<3; r++) {
268: valsrand[r]=5*x[r]*(1-x[r])*y[r]*(1-y[r]);
269: if (x[r]>=0.5 && y[r]>=0.5) {
270: vals1[r]=0.75;
271: vals2[r]=0.0;
272: }
273: if (x[r]>=0.5 && y[r]<0.5) {
274: vals1[r]=0.0;
275: vals2[r]=0.0;
276: }
277: if (x[r]<0.5 && y[r]>=0.5) {
278: vals1[r]=0.0;
279: vals2[r]=0.75;
280: }
281: if (x[r]<0.5 && y[r]<0.5) {
282: vals1[r]=0.75;
283: vals2[r]=0.0;
284: }
285: }
287: VecSetValues(user->work1,3,idx,vals1,INSERT_VALUES);
288: VecSetValues(user->work2,3,idx,vals2,INSERT_VALUES);
289: VecSetValues(user->work3,3,idx,valsrand,INSERT_VALUES);
290: }
292: VecAssemblyBegin(user->work1);
293: VecAssemblyEnd(user->work1);
294: VecAssemblyBegin(user->work2);
295: VecAssemblyEnd(user->work2);
296: VecAssemblyBegin(user->work3);
297: VecAssemblyEnd(user->work3);
299: VecAXPY(user->work1,1.0,user->work3);
300: VecAXPY(user->work2,1.0,user->work3);
302: for (i=0; i<n/4; i++) {
303: xx[4*i] = w1[i];
304: if (xx[4*i]>1) xx[4*i]=1;
306: xx[4*i+1] = w2[i];
307: if (xx[4*i+1]>1) xx[4*i+1]=1;
309: if (xx[4*i]+xx[4*i+1]>1) xx[4*i+1] = 1.0 - xx[4*i];
311: xx[4*i+2] = 1.0 - xx[4*i] - xx[4*i+1];
312: xx[4*i+3] = 0.0;
314: u1[i] = xx[4*i];
315: u2[i] = xx[4*i+1];
316: u3[i] = xx[4*i+2];
317: }
319: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
320: VecView(user->u1,view_out);
321: VecView(user->u2,view_out);
322: VecView(user->u3,view_out);
323: PetscViewerDestroy(&view_out);
325: DMDARestoreElements(user->da,&nele,&nen,&ele);
326: VecRestoreArrayRead(coords,&_coords);
327: VecRestoreArray(X,&xx);
328: VecRestoreArray(user->work2,&w1);
329: VecRestoreArray(user->work4,&w2);
330: VecRestoreArray(user->u1,&u1);
331: VecRestoreArray(user->u2,&u2);
332: VecRestoreArray(user->u3,&u3);
333: VecDestroy(&rand1);
334: VecDestroy(&rand2);
335: return(0);
336: }
342: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
343: {
345: AppCtx *user=(AppCtx*)ctx;
348: MatMultAdd(user->M,X,user->q,F);
349: return(0);
350: }
354: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat B,void *ctx)
355: {
357: AppCtx *user=(AppCtx*)ctx;
360: MatCopy(user->M,J,SAME_NONZERO_PATTERN);
361: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
362: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
363: return(0);
364: }
368: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
369: {
371: PetscScalar ***l,***u;
372: PetscInt xs,xm,ys,ym;
373: PetscInt j,i;
376: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
377: DMDAVecGetArrayDOF(da,xl,&l);
378: DMDAVecGetArrayDOF(da,xu,&u);
379: for (j=ys; j < ys+ym; j++) {
380: for (i=xs; i < xs+xm; i++) {
381: l[j][i][0] = 0.0;
382: l[j][i][1] = 0.0;
383: l[j][i][2] = 0.0;
384: l[j][i][3] = -PETSC_INFINITY;
385: u[j][i][0] = 1.0;
386: u[j][i][1] = 1.0;
387: u[j][i][2] = 1.0;
388: u[j][i][3] = PETSC_INFINITY;
389: }
390: }
391: DMDAVecRestoreArrayDOF(da,xl,&l);
392: DMDAVecRestoreArrayDOF(da,xu,&u);
393: return(0);
394: }
398: PetscErrorCode GetParams(AppCtx *user)
399: {
401: PetscBool flg;
404: /* Set default parameters */
405: user->xmin = 0.0; user->xmax = 1.0;
406: user->ymin = 0.0; user->ymax = 1.0;
407: user->T = 0.2; user->dt = 0.001;
408: user->epsilon = 0.05;
410: PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
411: PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
412: PetscOptionsGetReal(NULL,"-ymin",&user->ymin,&flg);
413: PetscOptionsGetReal(NULL,"-ymax",&user->ymax,&flg);
414: PetscOptionsGetReal(NULL,"-T",&user->T,&flg);
415: PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);
416: PetscOptionsGetScalar(NULL,"-epsilon",&user->epsilon,&flg);
417: return(0);
418: }
422: PetscErrorCode SetUpMatrices(AppCtx *user)
423: {
424: PetscErrorCode ierr;
425: PetscInt nele,nen,i,j;
426: const PetscInt *ele;
427: PetscScalar dt=user->dt;
428: PetscScalar y[3];
429: PetscInt idx[3];
430: PetscScalar eM_0[3][3],eM_2_odd[3][3],eM_2_even[3][3];
431: Mat M =user->M;
432: PetscScalar epsilon=user->epsilon;
433: PetscScalar hx;
434: PetscInt n,Mda,Nda;
435: DM da;
438: /* Create the mass matrix M_0 */
439: MatGetLocalSize(M,&n,NULL);
442: /* MatCreate(PETSC_COMM_WORLD,&user->M_0);*/
443: DMDAGetInfo(user->da,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
444: hx = 1.0/(Mda-1);
445: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,Mda,Nda,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
446: DMSetMatType(da,MATAIJ);
447: DMCreateMatrix(da,&user->M_0);
448: DMDestroy(&da);
450: eM_0[0][0]=eM_0[1][1]=eM_0[2][2]=hx*hx/12.0;
451: eM_0[0][1]=eM_0[0][2]=eM_0[1][0]=eM_0[1][2]=eM_0[2][0]=eM_0[2][1]=hx*hx/24.0;
453: eM_2_odd[0][0] = eM_2_odd[0][1] = eM_2_odd[0][2] = 0.0;
454: eM_2_odd[1][0] = eM_2_odd[1][1] = eM_2_odd[1][2] = 0.0;
455: eM_2_odd[2][0] = eM_2_odd[2][1] = eM_2_odd[2][2] = 0.0;
457: eM_2_odd[0][0]=1.0;
458: eM_2_odd[1][1]=eM_2_odd[2][2]=0.5;
459: eM_2_odd[0][1]=eM_2_odd[0][2]=eM_2_odd[1][0]=eM_2_odd[2][0]=-0.5;
461: eM_2_even[0][0] = eM_2_even[0][1] = eM_2_even[0][2] = 0.0;
462: eM_2_even[0][0] = eM_2_even[0][1] = eM_2_even[0][2] = 0.0;
463: eM_2_even[0][0] = eM_2_even[0][1] = eM_2_even[0][2] = 0.0;
465: eM_2_even[1][1]=1;
466: eM_2_even[0][0]=eM_2_even[2][2]=0.5;
467: eM_2_even[0][1]=eM_2_even[1][0]=eM_2_even[1][2]=eM_2_even[2][1]=-0.5;
469: /* Get local element info */
470: DMDAGetElements(user->da,&nele,&nen,&ele);
471: for (i=0; i < nele; i++) {
472: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
474: PetscInt row,cols[3],r,row_M_0;
475: PetscScalar vals[3],vals_M_0[3];
477: for (r=0; r<3; r++) {
478: row_M_0 = idx[r];
480: vals_M_0[0]=eM_0[r][0];
481: vals_M_0[1]=eM_0[r][1];
482: vals_M_0[2]=eM_0[r][2];
484: MatSetValues(user->M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);
486: if (y[1]==y[0]) {
487: row = 4*idx[r];
488: cols[0] = 4*idx[0]; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_odd[r][0];
489: cols[1] = 4*idx[1]; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_odd[r][1];
490: cols[2] = 4*idx[2]; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_odd[r][2];
491: /* Insert values in matrix M for 1st dof */
492: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
494: row = 4*idx[r]+1;
495: cols[0] = 4*idx[0]+1; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_odd[r][0];
496: cols[1] = 4*idx[1]+1; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_odd[r][1];
497: cols[2] = 4*idx[2]+1; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_odd[r][2];
498: /* Insert values in matrix M for 2nd dof */
499: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
501: row = 4*idx[r]+2;
502: cols[0] = 4*idx[0]+2; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_odd[r][0];
503: cols[1] = 4*idx[1]+2; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_odd[r][1];
504: cols[2] = 4*idx[2]+2; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_odd[r][2];
505: /* Insert values in matrix M for 3nd dof */
506: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
507: } else {
508: row = 4*idx[r];
509: cols[0] = 4*idx[0]; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_even[r][0];
510: cols[1] = 4*idx[1]; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_even[r][1];
511: cols[2] = 4*idx[2]; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_even[r][2];
512: /* Insert values in matrix M for 1st dof */
513: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
515: row = 4*idx[r]+1;
516: cols[0] = 4*idx[0]+1; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_even[r][0];
517: cols[1] = 4*idx[1]+1; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_even[r][1];
518: cols[2] = 4*idx[2]+1; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_even[r][2];
519: /* Insert values in matrix M for 2nd dof */
520: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
522: row = 4*idx[r]+2;
523: cols[0] = 4*idx[0]+2; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_even[r][0];
524: cols[1] = 4*idx[1]+2; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_even[r][1];
525: cols[2] = 4*idx[2]+2; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_even[r][2];
526: /* Insert values in matrix M for 3nd dof */
527: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
528: }
529: }
530: }
532: MatAssemblyBegin(user->M_0,MAT_FINAL_ASSEMBLY);
533: MatAssemblyEnd(user->M_0,MAT_FINAL_ASSEMBLY);
535: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
536: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
538: PetscScalar vals[9];
540: vals[0] = -1.0; vals[1] = 0.0; vals[2] = 0.0;
541: vals[3] = 0.0; vals[4] = -1.0; vals[5] = 0.0;
542: vals[6] = 0.0; vals[7] = 0.0; vals[8] = -1.0;
545: for (j=0; j < nele; j++) {
546: idx[0] = ele[3*j]; idx[1] = ele[3*j+1]; idx[2] = ele[3*j+2];
548: PetscInt r,rows[3],cols[3];
549: for (r=0; r<3; r++) {
551: rows[0] = 4*idx[0]+r; cols[0] = 4*idx[0]+3;
552: rows[1] = 4*idx[1]+r; cols[1] = 4*idx[1]+3;
553: rows[2] = 4*idx[2]+r; cols[2] = 4*idx[2]+3;
555: MatSetValuesLocal(M,3,rows,3,cols,vals,INSERT_VALUES);
556: MatSetValuesLocal(M,3,cols,3,rows,vals,INSERT_VALUES);
558: }
560: }
562: DMDARestoreElements(user->da,&nele,&nen,&ele);
564: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
565: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
569: VecCreate(PETSC_COMM_WORLD,&user->u1);
570: VecSetSizes(user->u1,n/4,PETSC_DECIDE);
571: VecSetFromOptions(user->u1);
572: VecDuplicate(user->u1,&user->u2);
573: VecDuplicate(user->u1,&user->u3);
574: VecDuplicate(user->u1,&user->work1);
575: VecDuplicate(user->u1,&user->work2);
576: VecDuplicate(user->u1,&user->work3);
577: VecDuplicate(user->u1,&user->work4);
578: return(0);
579: }