Actual source code: ex55.c

petsc-3.5.1 2014-07-24
Report Typos and Errors
  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: }