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: }