Actual source code: ex55.c
petsc-3.3-p7 2013-05-11
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)
17:
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*,MatStructure*,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);
58:
59: /* Get physics and time parameters */
60: GetParams(&user);
61: /* Create a 2D DA with dof = 2 */
62: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,PETSC_NULL,PETSC_NULL,&user.da);
63: /* Set Element type (triangular) */
64: DMDASetElementType(user.da,DMDA_ELEMENT_P1);
65:
66: /* Set x and y coordinates */
67: DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0.0,1.0);
68:
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);
75:
76: /* Get Jacobian matrix structure from the da */
77: DMCreateMatrix(user.da,MATAIJ,&user.M);
78: /* Form the jacobian matrix and M_0 */
79: SetUpMatrices(&user);
80: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
81:
82: /* Create nonlinear solver context */
83: SNESCreate(PETSC_COMM_WORLD,&snes);
84: SNESSetDM(snes,user.da);
85:
86: /* Set Function evaluation and jacobian evaluation routines */
87: SNESSetFunction(snes,r,FormFunction,(void*)&user);
88: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
89:
90: /* Set the boundary conditions */
91: SetVariableBounds(user.da,xl,xu);
92: SNESVISetVariableBounds(snes,xl,xu);
93: SNESSetFromOptions(snes);
94:
95: SetInitialGuess(x,&user);
96: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
97: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
98: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file1",FILE_MODE_WRITE,&view1);
99: /* Begin time loop */
100: while(t < user.T) {
101: VecView(user.u1,view1);
102: VecView(user.u2,view1);
103: VecView(user.u3,view1);
105: Update_q(user.q,user.u1,user.u2,user.u3,user.M_0,&user);
106: VecView(user.q,view_q);
107: SNESSolve(snes,PETSC_NULL,x);
108: PetscInt its;
109: SNESGetIterationNumber(snes,&its);
110: PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4f in %d iterations\n",t,its);
111: Update_u(user.u1,user.u2,user.u3,x);
112: t = t + user.dt;
113: VecView(user.u1,view_out);
114: VecView(user.u2,view_out);
115: VecView(user.u3,view_out);
116: }
117:
118:
119:
120: PetscViewerDestroy(&view_out);
121: PetscViewerDestroy(&view_q);
122: PetscViewerDestroy(&view1);
123:
124: VecDestroy(&x);
125: VecDestroy(&r);
126: VecDestroy(&xl);
127: VecDestroy(&xu);
128: VecDestroy(&user.q);
129: VecDestroy(&user.u1);
130: VecDestroy(&user.u2);
131: VecDestroy(&user.u3);
132: VecDestroy(&user.work1);
133: VecDestroy(&user.work2);
134: VecDestroy(&user.work3);
135: VecDestroy(&user.work4);
136: MatDestroy(&user.M);
137: MatDestroy(&user.M_0);
138: MatDestroy(&J);
139: DMDestroy(&user.da);
140: SNESDestroy(&snes);
141: PetscFinalize();
142: return 0;
143: }
147: PetscErrorCode Update_u(Vec u1,Vec u2,Vec u3,Vec X)
148: {
150: PetscInt i,n;
151: PetscScalar *u1_arr,*u2_arr,*u3_arr,*x_arr;
152:
154: VecGetLocalSize(u1,&n);
155: VecGetArray(u1,&u1_arr);
156: VecGetArray(u2,&u2_arr);
157: VecGetArray(u3,&u3_arr);
158: VecGetArray(X,&x_arr);
159: for(i=0;i<n;i++) {
160: u1_arr[i] = x_arr[4*i];
161: u2_arr[i] = x_arr[4*i+1];
162: u3_arr[i] = x_arr[4*i+2];
163: }
164: VecRestoreArray(u1,&u1_arr);
165: VecRestoreArray(u2,&u2_arr);
166: VecRestoreArray(u3,&u3_arr);
167: VecRestoreArray(X,&x_arr);
168: return(0);
169: }
173: PetscErrorCode Update_q(Vec q,Vec u1,Vec u2,Vec u3,Mat M_0,AppCtx *user)
174: {
176: PetscScalar *q_arr,*w_arr;
177: PetscInt i,n;
178: //PetscViewer view_q;
179:
181: VecSet(user->work1,user->dt/3);
182: // VecView(user->work1,PETSC_VIEWER_STDOUT_WORLD);
183: MatMult(M_0,user->work1,user->work2);
184: // VecView(user->work2,PETSC_VIEWER_STDOUT_WORLD);
185:
186: MatMult(M_0,u1,user->work1);
187: MatMult(M_0,u1,user->work4);
188: // VecView(u1,PETSC_VIEWER_STDOUT_WORLD);
189: // VecView(user->work4,PETSC_VIEWER_STDOUT_WORLD);
190: VecScale(user->work1,-1.0-(user->dt));
191: VecAXPY(user->work1,1.0,user->work2);
193: VecGetLocalSize(u1,&n);
194: VecGetArray(q,&q_arr);
195: VecGetArray(user->work1,&w_arr);
196: for(i=0;i<n;i++) {
197: q_arr[4*i] = w_arr[i];
198: }
199: VecRestoreArray(q,&q_arr);
200: VecRestoreArray(user->work1,&w_arr);
202: MatMult(M_0,u2,user->work1);
203: VecScale(user->work1,-1.0-(user->dt));
204: VecAXPY(user->work1,1.0,user->work2);
205:
206: VecGetArray(q,&q_arr);
207: VecGetArray(user->work1,&w_arr);
208: for(i=0;i<n;i++) {
209: q_arr[4*i+1] = w_arr[i];
210: }
211: VecRestoreArray(q,&q_arr);
212: VecRestoreArray(user->work1,&w_arr);
213:
214: MatMult(M_0,u3,user->work1);
215: VecScale(user->work1,-1.0-(user->dt));
216: VecAXPY(user->work1,1.0,user->work2);
218: VecGetArray(q,&q_arr);
219: VecGetArray(user->work1,&w_arr);
220: for(i=0;i<n;i++) {
221: q_arr[4*i+2] = w_arr[i];
222: q_arr[4*i+3] = 1.0;
223: }
224: VecRestoreArray(q,&q_arr);
225: VecRestoreArray(user->work1,&w_arr);
226: return(0);
227: }
231: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
232: {
233: PetscErrorCode ierr;
234: PetscInt nele,nen,n,i;
235: const PetscInt *ele;
236: Vec coords, rand1, rand2;
237: const PetscScalar *_coords;
238: PetscScalar x[3],y[3];
239: PetscInt idx[3];
240: PetscScalar *xx,*w1,*w2,*u1,*u2,*u3;
241: PetscViewer view_out;
244: /* Get ghosted coordinates */
245: DMDAGetGhostedCoordinates(user->da,&coords);
246: VecDuplicate(user->u1,&rand1);
247: VecDuplicate(user->u1,&rand2);
248: VecSetRandom(rand1,PETSC_NULL);
249: VecSetRandom(rand2,PETSC_NULL);
250:
251: VecGetLocalSize(X,&n);
252: VecGetArrayRead(coords,&_coords);
253: VecGetArray(X,&xx);
254: VecGetArray(user->work1,&w1);
255: VecGetArray(user->work2,&w2);
256: VecGetArray(user->u1,&u1);
257: VecGetArray(user->u2,&u2);
258: VecGetArray(user->u3,&u3);
259:
260: /* Get local element info */
261: DMDAGetElements(user->da,&nele,&nen,&ele);
262: for(i=0;i < nele;i++) {
263: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
264: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
265: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
266: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
267:
268: PetscScalar vals1[3],vals2[3],valsrand[3];
269: PetscInt r;
270: for(r=0;r<3;r++) {
271: valsrand[r]=5*x[r]*(1-x[r])*y[r]*(1-y[r]);
272: if (x[r]>=0.5 && y[r]>=0.5){
273: vals1[r]=0.75;
274: vals2[r]=0.0;
275: }
276: if (x[r]>=0.5 && y[r]<0.5){
277: vals1[r]=0.0;
278: vals2[r]=0.0;
279: }
280: if (x[r]<0.5 && y[r]>=0.5){
281: vals1[r]=0.0;
282: vals2[r]=0.75;
283: }
284: if (x[r]<0.5 && y[r]<0.5){
285: vals1[r]=0.75;
286: vals2[r]=0.0;
287: }
288: }
289:
290: VecSetValues(user->work1,3,idx,vals1,INSERT_VALUES);
291: VecSetValues(user->work2,3,idx,vals2,INSERT_VALUES);
292: VecSetValues(user->work3,3,idx,valsrand,INSERT_VALUES);
293: }
294:
295: VecAssemblyBegin(user->work1);
296: VecAssemblyEnd(user->work1);
297: VecAssemblyBegin(user->work2);
298: VecAssemblyEnd(user->work2);
299: VecAssemblyBegin(user->work3);
300: VecAssemblyEnd(user->work3);
301:
302: VecAXPY(user->work1,1.0,user->work3);
303: VecAXPY(user->work2,1.0,user->work3);
304:
305: for (i=0;i<n/4;i++) {
306: xx[4*i] = w1[i];
307: if (xx[4*i]>1) {
308: xx[4*i]=1;
309: }
310: xx[4*i+1] = w2[i];
311: if (xx[4*i+1]>1) {
312: xx[4*i+1]=1;
313: }
314: if (xx[4*i]+xx[4*i+1]>1){
315: xx[4*i+1] = 1.0 - xx[4*i];
316: }
317: xx[4*i+2] = 1.0 - xx[4*i] - xx[4*i+1];
318: xx[4*i+3] = 0.0;
319:
320: u1[i] = xx[4*i];
321: u2[i] = xx[4*i+1];
322: u3[i] = xx[4*i+2];
323: }
324:
325: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
326: VecView(user->u1,view_out);
327: VecView(user->u2,view_out);
328: VecView(user->u3,view_out);
329: PetscViewerDestroy(&view_out);
330:
331: DMDARestoreElements(user->da,&nele,&nen,&ele);
332: VecRestoreArrayRead(coords,&_coords);
333: VecRestoreArray(X,&xx);
334: VecRestoreArray(user->work2,&w1);
335: VecRestoreArray(user->work4,&w2);
336: VecRestoreArray(user->u1,&u1);
337: VecRestoreArray(user->u2,&u2);
338: VecRestoreArray(user->u3,&u3);
339: VecDestroy(&rand1);
340: VecDestroy(&rand2);
341: return(0);
342: }
348: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
349: {
351: AppCtx *user=(AppCtx*)ctx;
352:
354: MatMultAdd(user->M,X,user->q,F);
355: return(0);
356: }
360: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
361: {
363: AppCtx *user=(AppCtx*)ctx;
364:
366: *flg = SAME_NONZERO_PATTERN;
367: MatCopy(user->M,*J,*flg);
368: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
369: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
370: return(0);
371: }
375: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
376: {
378: PetscScalar ***l,***u;
379: PetscInt xs,xm,ys,ym;
380: PetscInt j,i;
381:
383:
384: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
385: DMDAVecGetArrayDOF(da,xl,&l);
386: DMDAVecGetArrayDOF(da,xu,&u);
387: for(j=ys; j < ys+ym; j++) {
388: for(i=xs; i < xs+xm;i++) {
389: l[j][i][0] = 0.0;
390: l[j][i][1] = 0.0;
391: l[j][i][2] = 0.0;
392: l[j][i][3] = -SNES_VI_INF;
393: u[j][i][0] = 1.0;
394: u[j][i][1] = 1.0;
395: u[j][i][2] = 1.0;
396: u[j][i][3] = SNES_VI_INF;
397: }
398: }
399: DMDAVecRestoreArrayDOF(da,xl,&l);
400: DMDAVecRestoreArrayDOF(da,xu,&u);
401: return(0);
402: }
406: PetscErrorCode GetParams(AppCtx* user)
407: {
409: PetscBool flg;
410:
412:
413: /* Set default parameters */
414: user->xmin = 0.0; user->xmax = 1.0;
415: user->ymin = 0.0; user->ymax = 1.0;
416: user->T = 0.2; user->dt = 0.001;
417: user->epsilon = 0.05;
418:
419: PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);
420: PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);
421: PetscOptionsGetReal(PETSC_NULL,"-ymin",&user->ymin,&flg);
422: PetscOptionsGetReal(PETSC_NULL,"-ymax",&user->ymax,&flg);
423: PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);
424: PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);
425: PetscOptionsGetScalar(PETSC_NULL,"-epsilon",&user->epsilon,&flg);
426: return(0);
427: }
431: PetscErrorCode SetUpMatrices(AppCtx* user)
432: {
433: PetscErrorCode ierr;
434: PetscInt nele,nen,i,j;
435: const PetscInt *ele;
436: PetscScalar dt=user->dt;
437: Vec coords;
438: const PetscScalar *_coords;
439: PetscScalar x[3],y[3];
440: PetscInt idx[3];
441: PetscScalar eM_0[3][3],eM_2_odd[3][3],eM_2_even[3][3];
442: Mat M=user->M;
443: PetscScalar epsilon=user->epsilon;
444: PetscScalar hx;
445: PetscInt n,Mda,Nda;
446: DM da;
447:
449: /* Get ghosted coordinates */
450: DMDAGetGhostedCoordinates(user->da,&coords);
451: VecGetArrayRead(coords,&_coords);
452:
453: /* Create the mass matrix M_0 */
454: MatGetLocalSize(M,&n,PETSC_NULL);
455:
456:
457: /* MatCreate(PETSC_COMM_WORLD,&user->M_0);*/
458: DMDAGetInfo(user->da,PETSC_NULL,&Mda,&Nda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
459: hx = 1.0/(Mda-1);
460: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,Mda,Nda,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da);
461: DMCreateMatrix(da,MATAIJ,&user->M_0);
462: DMDestroy(&da);
463:
464: eM_0[0][0]=eM_0[1][1]=eM_0[2][2]=hx*hx/12.0;
465: 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;
466:
467: eM_2_odd[0][0] = eM_2_odd[0][1] = eM_2_odd[0][2] = 0.0;
468: eM_2_odd[1][0] = eM_2_odd[1][1] = eM_2_odd[1][2] = 0.0;
469: eM_2_odd[2][0] = eM_2_odd[2][1] = eM_2_odd[2][2] = 0.0;
471: eM_2_odd[0][0]=1.0;
472: eM_2_odd[1][1]=eM_2_odd[2][2]=0.5;
473: eM_2_odd[0][1]=eM_2_odd[0][2]=eM_2_odd[1][0]=eM_2_odd[2][0]=-0.5;
474:
475: eM_2_even[0][0] = eM_2_even[0][1] = eM_2_even[0][2] = 0.0;
476: eM_2_even[0][0] = eM_2_even[0][1] = eM_2_even[0][2] = 0.0;
477: eM_2_even[0][0] = eM_2_even[0][1] = eM_2_even[0][2] = 0.0;
478:
479: eM_2_even[1][1]=1;
480: eM_2_even[0][0]=eM_2_even[2][2]=0.5;
481: eM_2_even[0][1]=eM_2_even[1][0]=eM_2_even[1][2]=eM_2_even[2][1]=-0.5;
483: /* Get local element info */
484: DMDAGetElements(user->da,&nele,&nen,&ele);
485: for(i=0;i < nele;i++) {
486: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
487: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
488: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
489: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
490:
491: PetscInt row,cols[3],r,row_M_0;
492: PetscScalar vals[3],vals_M_0[3];
493:
494: for(r=0;r<3;r++) {
495: row_M_0 = idx[r];
496:
497: vals_M_0[0]=eM_0[r][0];
498: vals_M_0[1]=eM_0[r][1];
499: vals_M_0[2]=eM_0[r][2];
500:
501: MatSetValues(user->M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);
502:
503: if (y[1]==y[0]) {
504: row = 4*idx[r];
505: cols[0] = 4*idx[0]; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_odd[r][0];
506: cols[1] = 4*idx[1]; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_odd[r][1];
507: cols[2] = 4*idx[2]; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_odd[r][2];
508: /* Insert values in matrix M for 1st dof */
509: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
510:
511: row = 4*idx[r]+1;
512: cols[0] = 4*idx[0]+1; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_odd[r][0];
513: cols[1] = 4*idx[1]+1; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_odd[r][1];
514: cols[2] = 4*idx[2]+1; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_odd[r][2];
515: /* Insert values in matrix M for 2nd dof */
516: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
517:
518: row = 4*idx[r]+2;
519: cols[0] = 4*idx[0]+2; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_odd[r][0];
520: cols[1] = 4*idx[1]+2; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_odd[r][1];
521: cols[2] = 4*idx[2]+2; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_odd[r][2];
522: /* Insert values in matrix M for 3nd dof */
523: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
524: }
525: else{
526: row = 4*idx[r];
527: cols[0] = 4*idx[0]; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_even[r][0];
528: cols[1] = 4*idx[1]; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_even[r][1];
529: cols[2] = 4*idx[2]; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_even[r][2];
530: /* Insert values in matrix M for 1st dof */
531: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
532:
533: row = 4*idx[r]+1;
534: cols[0] = 4*idx[0]+1; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_even[r][0];
535: cols[1] = 4*idx[1]+1; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_even[r][1];
536: cols[2] = 4*idx[2]+1; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_even[r][2];
537: /* Insert values in matrix M for 2nd dof */
538: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
539:
540: row = 4*idx[r]+2;
541: cols[0] = 4*idx[0]+2; vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_even[r][0];
542: cols[1] = 4*idx[1]+2; vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_even[r][1];
543: cols[2] = 4*idx[2]+2; vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_even[r][2];
544: /* Insert values in matrix M for 3nd dof */
545: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
546: }
547: }
548: }
549:
550: MatAssemblyBegin(user->M_0,MAT_FINAL_ASSEMBLY);
551: MatAssemblyEnd(user->M_0,MAT_FINAL_ASSEMBLY);
552:
553: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
554: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
555:
556: PetscScalar vals[9];
557:
558: vals[0] = -1.0; vals[1] = 0.0; vals[2] = 0.0;
559: vals[3] = 0.0; vals[4] = -1.0; vals[5] = 0.0;
560: vals[6] = 0.0; vals[7] = 0.0; vals[8] = -1.0;
561:
562:
563: for(j=0;j < nele;j++) {
564: idx[0] = ele[3*j]; idx[1] = ele[3*j+1]; idx[2] = ele[3*j+2];
565:
566: PetscInt r,rows[3],cols[3];
567: for(r=0;r<3;r++) {
568:
569: rows[0] = 4*idx[0]+r; cols[0] = 4*idx[0]+3;
570: rows[1] = 4*idx[1]+r; cols[1] = 4*idx[1]+3;
571: rows[2] = 4*idx[2]+r; cols[2] = 4*idx[2]+3;
572: MatSetValuesLocal(M,3,rows,3,cols,vals,INSERT_VALUES);
573: MatSetValuesLocal(M,3,cols,3,rows,vals,INSERT_VALUES);
574:
575: }
576:
577: }
578:
579: DMDARestoreElements(user->da,&nele,&nen,&ele);
580: VecRestoreArrayRead(coords,&_coords);
581:
582: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
583: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
584:
585:
586:
587: VecCreate(PETSC_COMM_WORLD,&user->u1);
588: VecSetSizes(user->u1,n/4,PETSC_DECIDE);
589: VecSetFromOptions(user->u1);
590: VecDuplicate(user->u1,&user->u2);
591: VecDuplicate(user->u1,&user->u3);
592: VecDuplicate(user->u1,&user->work1);
593: VecDuplicate(user->u1,&user->work2);
594: VecDuplicate(user->u1,&user->work3);
595: VecDuplicate(user->u1,&user->work4);
596:
597: return(0);
598: }