Actual source code: ex55.c

  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:         DMGetMatrix(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:         SNESSetType(snes,SNESVI);
 91:         SNESSetFromOptions(snes);
 92:         /* Set the boundary conditions */
 93:         SetVariableBounds(user.da,xl,xu);
 94:         SNESVISetVariableBounds(snes,xl,xu);
 95: 
 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,PETSC_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:         }
118: 
119: 
120: 
121:         PetscViewerDestroy(&view_out);
122:         PetscViewerDestroy(&view_q);
123:         PetscViewerDestroy(&view1);
124: 
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;
153: 
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;
180: 
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);
186: 
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++) {
198:     q_arr[4*i] = w_arr[i];
199:   }
200:   VecRestoreArray(q,&q_arr);
201:   VecRestoreArray(user->work1,&w_arr);

203:   MatMult(M_0,u2,user->work1);
204:   VecScale(user->work1,-1.0-(user->dt));
205:   VecAXPY(user->work1,1.0,user->work2);
206: 
207:   VecGetArray(q,&q_arr);
208:   VecGetArray(user->work1,&w_arr);
209:   for(i=0;i<n;i++) {
210:     q_arr[4*i+1] = w_arr[i];
211:   }
212:   VecRestoreArray(q,&q_arr);
213:   VecRestoreArray(user->work1,&w_arr);
214: 
215:   MatMult(M_0,u3,user->work1);
216:   VecScale(user->work1,-1.0-(user->dt));
217:   VecAXPY(user->work1,1.0,user->work2);

219:   VecGetArray(q,&q_arr);
220:   VecGetArray(user->work1,&w_arr);
221:     for(i=0;i<n;i++) {
222:     q_arr[4*i+2] = w_arr[i];
223:     q_arr[4*i+3] = 1.0;
224:   }
225:   VecRestoreArray(q,&q_arr);
226:   VecRestoreArray(user->work1,&w_arr);
227:   return(0);
228: }

232: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
233: {
234:         PetscErrorCode    ierr;
235:         PetscInt          nele,nen,n,i;
236:         const PetscInt    *ele;
237:         Vec               coords, rand1, rand2;
238:         const PetscScalar *_coords;
239:         PetscScalar       x[3],y[3];
240:         PetscInt          idx[3];
241:         PetscScalar           *xx,*w1,*w2,*u1,*u2,*u3;
242:         PetscViewer                  view_out;

245:         /* Get ghosted coordinates */
246:         DMDAGetGhostedCoordinates(user->da,&coords);
247:         VecDuplicate(user->u1,&rand1);
248:         VecDuplicate(user->u1,&rand2);
249:         VecSetRandom(rand1,PETSC_NULL);
250:         VecSetRandom(rand2,PETSC_NULL);
251: 
252:         VecGetLocalSize(X,&n);
253:         VecGetArrayRead(coords,&_coords);
254:         VecGetArray(X,&xx);
255:         VecGetArray(user->work1,&w1);
256:         VecGetArray(user->work2,&w2);
257:         VecGetArray(user->u1,&u1);
258:         VecGetArray(user->u2,&u2);
259:         VecGetArray(user->u3,&u3);
260: 
261:         /* Get local element info */
262:         DMDAGetElements(user->da,&nele,&nen,&ele);
263:         for(i=0;i < nele;i++) {
264:                 idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
265:                 x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
266:                 x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
267:                 x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
268: 
269:                 PetscScalar vals1[3],vals2[3],valsrand[3];
270:                 PetscInt r;
271:                 for(r=0;r<3;r++) {
272:                         valsrand[r]=5*x[r]*(1-x[r])*y[r]*(1-y[r]);
273:                         if (x[r]>=0.5 && y[r]>=0.5){
274:                                 vals1[r]=0.75;
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.0;
280:                         }
281:                         if (x[r]<0.5 && y[r]>=0.5){
282:                                 vals1[r]=0.0;
283:                                 vals2[r]=0.75;
284:                         }
285:                         if (x[r]<0.5 && y[r]<0.5){
286:                                 vals1[r]=0.75;
287:                                 vals2[r]=0.0;
288:                         }
289:                 }
290: 
291:                 VecSetValues(user->work1,3,idx,vals1,INSERT_VALUES);
292:                 VecSetValues(user->work2,3,idx,vals2,INSERT_VALUES);
293:                 VecSetValues(user->work3,3,idx,valsrand,INSERT_VALUES);
294:         }
295: 
296:         VecAssemblyBegin(user->work1);
297:         VecAssemblyEnd(user->work1);
298:         VecAssemblyBegin(user->work2);
299:         VecAssemblyEnd(user->work2);
300:         VecAssemblyBegin(user->work3);
301:         VecAssemblyEnd(user->work3);
302: 
303:         VecAXPY(user->work1,1.0,user->work3);
304:         VecAXPY(user->work2,1.0,user->work3);
305: 
306:         for (i=0;i<n/4;i++) {
307:                 xx[4*i] = w1[i];
308:                 if (xx[4*i]>1) {
309:                         xx[4*i]=1;
310:                 }
311:                 xx[4*i+1] = w2[i];
312:                 if (xx[4*i+1]>1) {
313:                         xx[4*i+1]=1;
314:                 }
315:                 if (xx[4*i]+xx[4*i+1]>1){
316:                         xx[4*i+1] = 1.0 - xx[4*i];
317:                 }
318:                 xx[4*i+2] = 1.0 - xx[4*i] - xx[4*i+1];
319:                 xx[4*i+3] = 0.0;
320: 
321:                 u1[i] = xx[4*i];
322:                 u2[i] = xx[4*i+1];
323:                 u3[i] = xx[4*i+2];
324:         }
325: 
326:         PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
327:         VecView(user->u1,view_out);
328:         VecView(user->u2,view_out);
329:         VecView(user->u3,view_out);
330:         PetscViewerDestroy(&view_out);
331: 
332:         DMDARestoreElements(user->da,&nele,&nen,&ele);
333:         VecRestoreArrayRead(coords,&_coords);
334:         VecRestoreArray(X,&xx);
335:         VecRestoreArray(user->work2,&w1);
336:         VecRestoreArray(user->work4,&w2);
337:         VecRestoreArray(user->u1,&u1);
338:         VecRestoreArray(user->u2,&u2);
339:         VecRestoreArray(user->u3,&u3);
340:         VecDestroy(&rand1);
341:         VecDestroy(&rand2);
342:         return(0);
343: }



349: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
350: {
352:   AppCtx         *user=(AppCtx*)ctx;
353: 
355:   MatMultAdd(user->M,X,user->q,F);
356:   return(0);
357: }

361: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
362: {
364:   AppCtx         *user=(AppCtx*)ctx;
365: 
367:   *flg = SAME_NONZERO_PATTERN;
368:   MatCopy(user->M,*J,*flg);
369:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
370:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
371:   return(0);
372: }

376: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
377: {
379:   PetscScalar    ***l,***u;
380:   PetscInt       xs,xm,ys,ym;
381:   PetscInt       j,i;
382: 
384: 
385:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
386:   DMDAVecGetArrayDOF(da,xl,&l);
387:   DMDAVecGetArrayDOF(da,xu,&u);
388:   for(j=ys; j < ys+ym; j++) {
389:     for(i=xs; i < xs+xm;i++) {
390:       l[j][i][0] = 0.0;
391:       l[j][i][1] = 0.0;
392:       l[j][i][2] = 0.0;
393:       l[j][i][3] = -SNES_VI_INF;
394:       u[j][i][0] = 1.0;
395:       u[j][i][1] = 1.0;
396:       u[j][i][2] = 1.0;
397:       u[j][i][3] = SNES_VI_INF;
398:     }
399:   }
400:   DMDAVecRestoreArrayDOF(da,xl,&l);
401:   DMDAVecRestoreArrayDOF(da,xu,&u);
402:   return(0);
403: }

407: PetscErrorCode GetParams(AppCtx* user)
408: {
410:   PetscBool      flg;
411: 
413: 
414:   /* Set default parameters */
415:   user->xmin = 0.0; user->xmax = 1.0;
416:   user->ymin = 0.0; user->ymax = 1.0;
417:   user->T = 0.2;    user->dt = 0.001;
418:   user->epsilon = 0.05;
419: 
420:   PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);
421:   PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);
422:   PetscOptionsGetReal(PETSC_NULL,"-ymin",&user->ymin,&flg);
423:   PetscOptionsGetReal(PETSC_NULL,"-ymax",&user->ymax,&flg);
424:   PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);
425:   PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);
426:   PetscOptionsGetScalar(PETSC_NULL,"-epsilon",&user->epsilon,&flg);
427:   return(0);
428: }

432: PetscErrorCode SetUpMatrices(AppCtx* user)
433: {
434:   PetscErrorCode    ierr;
435:   PetscInt          nele,nen,i,j;
436:   const PetscInt    *ele;
437:   PetscScalar       dt=user->dt;
438:   Vec               coords;
439:   const PetscScalar *_coords;
440:   PetscScalar       x[3],y[3];
441:   PetscInt          idx[3];
442:   PetscScalar       eM_0[3][3],eM_2_odd[3][3],eM_2_even[3][3];
443:   Mat               M=user->M;
444:   PetscScalar       epsilon=user->epsilon;
445:   PetscScalar                  hx;
446:   PetscInt n,Mda,Nda;
447:   DM               da;
448: 
450:   /* Get ghosted coordinates */
451:   DMDAGetGhostedCoordinates(user->da,&coords);
452:   VecGetArrayRead(coords,&_coords);
453: 
454:   /* Create the mass matrix M_0 */
455:   MatGetLocalSize(M,&n,PETSC_NULL);
456: 
457: 
458:   /* MatCreate(PETSC_COMM_WORLD,&user->M_0);*/
459:   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);
460:   hx = 1.0/(Mda-1);
461:   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);
462:   DMGetMatrix(da,MATAIJ,&user->M_0);
463:   DMDestroy(&da);
464: 
465:   eM_0[0][0]=eM_0[1][1]=eM_0[2][2]=hx*hx/12.0;
466:   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;
467: 
468:   eM_2_odd[0][0] = eM_2_odd[0][1] = eM_2_odd[0][2] = 0.0;
469:   eM_2_odd[1][0] = eM_2_odd[1][1] = eM_2_odd[1][2] = 0.0;
470:   eM_2_odd[2][0] = eM_2_odd[2][1] = eM_2_odd[2][2] = 0.0;

472:   eM_2_odd[0][0]=1.0;
473:   eM_2_odd[1][1]=eM_2_odd[2][2]=0.5;
474:   eM_2_odd[0][1]=eM_2_odd[0][2]=eM_2_odd[1][0]=eM_2_odd[2][0]=-0.5;
475: 
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:   eM_2_even[0][0] = eM_2_even[0][1] = eM_2_even[0][2] = 0.0;
479: 
480:   eM_2_even[1][1]=1;
481:   eM_2_even[0][0]=eM_2_even[2][2]=0.5;
482:   eM_2_even[0][1]=eM_2_even[1][0]=eM_2_even[1][2]=eM_2_even[2][1]=-0.5;

484:   /* Get local element info */
485:   DMDAGetElements(user->da,&nele,&nen,&ele);
486:   for(i=0;i < nele;i++) {
487:     idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
488:     x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
489:     x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
490:     x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
491: 
492:     PetscInt    row,cols[3],r,row_M_0;
493:     PetscScalar vals[3],vals_M_0[3];
494: 
495:     for(r=0;r<3;r++) {
496:       row_M_0 = idx[r];
497: 
498:       vals_M_0[0]=eM_0[r][0];
499:       vals_M_0[1]=eM_0[r][1];
500:       vals_M_0[2]=eM_0[r][2];
501: 
502:       MatSetValues(user->M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);
503: 
504:       if (y[1]==y[0]) {
505:         row = 4*idx[r];
506:         cols[0] = 4*idx[0];     vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_odd[r][0];
507:         cols[1] = 4*idx[1];     vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_odd[r][1];
508:         cols[2] = 4*idx[2];     vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_odd[r][2];
509:         /* Insert values in matrix M for 1st dof */
510:         MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
511: 
512:         row = 4*idx[r]+1;
513:         cols[0] = 4*idx[0]+1;   vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_odd[r][0];
514:         cols[1] = 4*idx[1]+1;   vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_odd[r][1];
515:         cols[2] = 4*idx[2]+1;   vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_odd[r][2];
516:         /* Insert values in matrix M for 2nd dof */
517:         MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
518: 
519:         row = 4*idx[r]+2;
520:         cols[0] = 4*idx[0]+2;   vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_odd[r][0];
521:         cols[1] = 4*idx[1]+2;   vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_odd[r][1];
522:         cols[2] = 4*idx[2]+2;   vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_odd[r][2];
523:         /* Insert values in matrix M for 3nd dof */
524:         MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
525:       }
526:       else{
527:         row = 4*idx[r];
528:         cols[0] = 4*idx[0];     vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_even[r][0];
529:         cols[1] = 4*idx[1];     vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_even[r][1];
530:         cols[2] = 4*idx[2];     vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_even[r][2];
531:         /* Insert values in matrix M for 1st dof */
532:         MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
533: 
534:         row = 4*idx[r]+1;
535:         cols[0] = 4*idx[0]+1;   vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_even[r][0];
536:         cols[1] = 4*idx[1]+1;   vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_even[r][1];
537:         cols[2] = 4*idx[2]+1;   vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_even[r][2];
538:         /* Insert values in matrix M for 2nd dof */
539:         MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
540: 
541:         row = 4*idx[r]+2;
542:         cols[0] = 4*idx[0]+2;   vals[0] = eM_0[r][0]+dt*epsilon*epsilon*eM_2_even[r][0];
543:         cols[1] = 4*idx[1]+2;   vals[1] = eM_0[r][1]+dt*epsilon*epsilon*eM_2_even[r][1];
544:         cols[2] = 4*idx[2]+2;   vals[2] = eM_0[r][2]+dt*epsilon*epsilon*eM_2_even[r][2];
545:         /* Insert values in matrix M for 3nd dof */
546:         MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
547:       }
548:     }
549:   }
550: 
551:   MatAssemblyBegin(user->M_0,MAT_FINAL_ASSEMBLY);
552:   MatAssemblyEnd(user->M_0,MAT_FINAL_ASSEMBLY);
553: 
554:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
555:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
556: 
557:   PetscScalar vals[9];
558: 
559:   vals[0] = -1.0; vals[1] =  0.0; vals[2] =  0.0;
560:   vals[3] =  0.0; vals[4] = -1.0; vals[5] =  0.0;
561:   vals[6] =  0.0; vals[7] =  0.0; vals[8] = -1.0;
562: 
563: 
564:   for(j=0;j < nele;j++) {
565:     idx[0] = ele[3*j]; idx[1] = ele[3*j+1]; idx[2] = ele[3*j+2];
566: 
567:     PetscInt   r,rows[3],cols[3];
568:     for(r=0;r<3;r++) {
569: 
570:       rows[0] = 4*idx[0]+r;        cols[0] = 4*idx[0]+3;
571:       rows[1] = 4*idx[1]+r;   cols[1] = 4*idx[1]+3;
572:       rows[2] = 4*idx[2]+r;   cols[2] = 4*idx[2]+3;
573:       MatSetValuesLocal(M,3,rows,3,cols,vals,INSERT_VALUES);
574:       MatSetValuesLocal(M,3,cols,3,rows,vals,INSERT_VALUES);
575: 
576:     }
577: 
578:   }
579: 
580:   DMDARestoreElements(user->da,&nele,&nen,&ele);
581:   VecRestoreArrayRead(coords,&_coords);
582: 
583:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
584:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
585: 
586: 
587: 
588:   VecCreate(PETSC_COMM_WORLD,&user->u1);
589:   VecSetSizes(user->u1,n/4,PETSC_DECIDE);
590:   VecSetFromOptions(user->u1);
591:   VecDuplicate(user->u1,&user->u2);
592:   VecDuplicate(user->u1,&user->u3);
593:   VecDuplicate(user->u1,&user->work1);
594:   VecDuplicate(user->u1,&user->work2);
595:   VecDuplicate(user->u1,&user->work3);
596:   VecDuplicate(user->u1,&user->work4);
597: 
598:   return(0);
599: }