Actual source code: ex64.c
petsc-3.4.2 2013-07-02

      1: static char help[] = "1D coupled Allen-Cahn and Cahn-Hilliard equation for constant mobility. Only c_v and eta are considered.\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: ./ex63 -ksp_type fgmres -snes_vi_monitor   -snes_atol 1.e-11 -snes_converged_reason -ksp_converged_reason   -snes_linesearch_monitor -VG 10000000 -draw_fields 1,3,4  -pc_type mg -pc_mg_galerkin -log_summary -dt .000001 -mg_coarse_pc_type svd  -ksp_monitor_true_residual -ksp_rtol 1.e-9 -snes_linesearch_type basic -T .0020 -P_casc .0005 -snes_monitor_solution -da_refine 10
     13:  */

     15:  #include petscsnes.h
     16:  #include petscdmda.h

     18: typedef struct {
     19:   PetscReal   dt,T; /* Time step and end time */
     20:   DM          da1,da1_clone,da2;
     21:   Mat         M;    /* Jacobian matrix */
     22:   Mat         M_0;
     23:   Vec         q,wv,cv,eta,DPsiv,DPsieta,logcv,logcv2;
     24:   Vec         work1,work2;
     25:   PetscScalar Mv,L,kaeta,kav,Evf,A,B,cv0,Sv;
     26:   PetscReal   xmin,xmax;
     27:   PetscInt    nx;
     28:   PetscBool   graphics;
     29:   PetscBool   periodic;
     30:   PetscBool   lumpedmass;
     31: } AppCtx;

     33: PetscErrorCode GetParams(AppCtx*);
     34: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
     35: PetscErrorCode SetUpMatrices(AppCtx*);
     36: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
     37: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
     38: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
     39: PetscErrorCode Update_q(AppCtx*);
     40: PetscErrorCode Update_u(Vec,AppCtx*);
     41: PetscErrorCode DPsi(AppCtx*);
     42: PetscErrorCode Llog(Vec,Vec);
     43: PetscErrorCode CheckRedundancy(SNES,IS,IS*,DM);

     47: int main(int argc, char **argv)
     48: {
     50:   Vec            x,r;  /* solution and residual vectors */
     51:   SNES           snes; /* Nonlinear solver context */
     52:   AppCtx         user; /* Application context */
     53:   Vec            xl,xu; /* Upper and lower bounds on variables */
     54:   Mat            J;
     55:   PetscScalar    t=0.0;
     56:   PetscViewer    view_out, view_p, view_q, view_psi, view_mat;
     57:   PetscReal      bounds[] = {1000.0,-1000.,0.0,1.0,1000.0,-1000.0,0.0,1.0,1000.0,-1000.0};

     60:   PetscInitialize(&argc,&argv, (char*)0, help);

     62:   /* Get physics and time parameters */
     63:   GetParams(&user);

     65:   if (user.periodic) {
     66:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 3, 1,NULL,&user.da1);
     67:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 3, 1,NULL,&user.da1_clone);
     68:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 1, 1,NULL,&user.da2);
     69:   } else {
     70:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE, -4, 3, 1,NULL,&user.da1);
     71:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE, -4, 3, 1,NULL,&user.da1_clone);
     72:     DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE, -4, 1, 1,NULL,&user.da2);

     74:   }
     75:   /* Set Element type (triangular) */
     76:   DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
     77:   DMDASetElementType(user.da2,DMDA_ELEMENT_P1);

     79:   /* Set x and y coordinates */
     80:   DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,NULL,NULL,NULL,NULL);
     81:   DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,NULL,NULL,NULL,NULL);
     82:   /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
     83:   DMCreateGlobalVector(user.da1,&x);
     84:   VecDuplicate(x,&r);
     85:   VecDuplicate(x,&xl);
     86:   VecDuplicate(x,&xu);
     87:   VecDuplicate(x,&user.q);

     89:   /* Get global vector user->wv from da2 and duplicate other vectors */
     90:   DMCreateGlobalVector(user.da2,&user.wv);
     91:   VecDuplicate(user.wv,&user.cv);
     92:   VecDuplicate(user.wv,&user.eta);
     93:   VecDuplicate(user.wv,&user.DPsiv);
     94:   VecDuplicate(user.wv,&user.DPsieta);
     95:   VecDuplicate(user.wv,&user.logcv);
     96:   VecDuplicate(user.wv,&user.logcv2);
     97:   VecDuplicate(user.wv,&user.work1);
     98:   VecDuplicate(user.wv,&user.work2);

    100:   /* Get Jacobian matrix structure from the da for the entire thing, da1 */
    101:   DMCreateMatrix(user.da1,MATAIJ,&user.M);
    102:   /* Get the (usual) mass matrix structure from da2 */
    103:   DMCreateMatrix(user.da2,MATAIJ,&user.M_0);
    104:   /* Form the jacobian matrix and M_0 */
    105:   SetUpMatrices(&user);
    106:   SetInitialGuess(x,&user);
    107:   MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);

    109:   SNESCreate(PETSC_COMM_WORLD,&snes);
    110:   SNESSetDM(snes,user.da1);

    112:   SNESSetFunction(snes,r,FormFunction,(void*)&user);
    113:   SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);

    115:   SetVariableBounds(user.da1,xl,xu);
    116:   SNESVISetVariableBounds(snes,xl,xu);
    117:   SNESSetFromOptions(snes);
    118:   /*SNESVISetRedundancyCheck(snes,(PetscErrorCode (*)(SNES,IS,IS*,void*))CheckRedundancy,user.da1_clone);*/
    119:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
    120:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_p",FILE_MODE_WRITE,&view_p);
    121:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
    122:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
    123:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
    124:   /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */

    126:   if (user.graphics) {
    127:     PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds);

    129:     VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
    130:   }
    131:   while (t<user.T) {

    133:     char        filename[PETSC_MAX_PATH_LEN];
    134:     PetscScalar a = 1.0;
    135:     PetscInt    i;
    136:     PetscViewer view;

    139:     SNESSetFunction(snes,r,FormFunction,(void*)&user);
    140:     SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);

    142:     DPsi(&user);
    143:     VecView(user.DPsiv,view_psi);
    144:     VecView(user.DPsieta,view_psi);

    146:     Update_q(&user);
    147:     VecView(user.q,view_q);
    148:     MatView(user.M,view_mat);

    150:     SNESSolve(snes,NULL,x);
    151:     VecView(x,view_out);

    154:     if (user.graphics) {
    155:       VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
    156:     }

    158:     PetscInt its;
    159:     SNESGetIterationNumber(snes,&its);
    160:     PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4g in %d iterations\n",t,its);

    162:     Update_u(x,&user);

    164:     for (i=0; i < (int)(user.T/a) ; i++) {
    165:       if (t/a > i - user.dt/a && t/a < i + user.dt/a) {
    166:         sprintf(filename,"output_%f",t);
    167:         PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&view);
    168:         VecView(user.cv,view);
    169:         VecView(user.eta,view);
    170:         PetscViewerDestroy(&view);
    171:       }
    172:     }

    174:     t = t + user.dt;

    176:   }

    179:   PetscViewerDestroy(&view_out);
    180:   PetscViewerDestroy(&view_p);
    181:   PetscViewerDestroy(&view_q);
    182:   PetscViewerDestroy(&view_psi);
    183:   PetscViewerDestroy(&view_mat);
    184:   VecDestroy(&x);
    185:   VecDestroy(&r);
    186:   VecDestroy(&xl);
    187:   VecDestroy(&xu);
    188:   VecDestroy(&user.q);
    189:   VecDestroy(&user.wv);
    190:   VecDestroy(&user.cv);
    191:   VecDestroy(&user.eta);
    192:   VecDestroy(&user.DPsiv);
    193:   VecDestroy(&user.DPsieta);
    194:   VecDestroy(&user.logcv);
    195:   VecDestroy(&user.logcv2);
    196:   VecDestroy(&user.work1);
    197:   VecDestroy(&user.work2);
    198:   MatDestroy(&user.M);
    199:   MatDestroy(&user.M_0);
    200:   DMDestroy(&user.da1);
    201:   DMDestroy(&user.da1_clone);
    202:   DMDestroy(&user.da2);
    203:   SNESDestroy(&snes);
    204:   PetscFinalize();
    205:   return 0;
    206: }

    210: PetscErrorCode Update_u(Vec X,AppCtx *user)
    211: {
    213:   PetscInt       i,n;
    214:   PetscScalar    *xx,*wv_p,*cv_p,*eta_p;

    217:   VecGetLocalSize(user->wv,&n);

    219:   VecGetArray(X,&xx);
    220:   VecGetArray(user->wv,&wv_p);
    221:   VecGetArray(user->cv,&cv_p);
    222:   VecGetArray(user->eta,&eta_p);

    225:   for (i=0; i<n; i++) {
    226:     wv_p[i]  = xx[3*i];
    227:     cv_p[i]  = xx[3*i+1];
    228:     eta_p[i] = xx[3*i+2];
    229:   }
    230:   VecRestoreArray(X,&xx);
    231:   VecRestoreArray(user->wv,&wv_p);
    232:   VecRestoreArray(user->cv,&cv_p);
    233:   VecRestoreArray(user->eta,&eta_p);
    234:   return(0);
    235: }

    239: PetscErrorCode Update_q(AppCtx *user)
    240: {
    242:   PetscScalar    *q_p, *w1, *w2;
    243:   PetscInt       n,i;

    246:   VecGetArray(user->q,&q_p);
    247:   VecGetArray(user->work1,&w1);
    248:   VecGetArray(user->work2,&w2);
    249:   VecGetLocalSize(user->wv,&n);

    251:   MatMult(user->M_0,user->cv,user->work1);
    252:   VecScale(user->work1,-1.0);
    253:   for (i=0; i<n; i++) q_p[3*i]=w1[i];

    255:   MatMult(user->M_0,user->DPsiv,user->work1);
    256:   for (i=0; i<n; i++) q_p[3*i+1]=w1[i];

    258:   VecCopy(user->DPsieta,user->work1);
    259:   VecScale(user->work1,user->L*user->dt);
    260:   VecAXPY(user->work1,-1.0,user->eta);
    261:   MatMult(user->M_0,user->work1,user->work2);
    262:   for (i=0; i<n; i++) q_p[3*i+2]=w2[i];

    264:   VecRestoreArray(user->q,&q_p);
    265:   VecRestoreArray(user->work1,&w1);
    266:   VecRestoreArray(user->work2,&w2);
    267:   return(0);
    268: }

    272: PetscErrorCode DPsi(AppCtx *user)
    273: {
    275:   PetscScalar    Evf=user->Evf,A=user->A,B=user->B,cv0=user->cv0;
    276:   PetscScalar    *cv_p,*eta_p,*logcv_p,*logcv2_p,*DPsiv_p,*DPsieta_p;
    277:   PetscInt       n,i;

    280:   VecGetLocalSize(user->cv,&n);
    281:   VecGetArray(user->cv,&cv_p);
    282:   VecGetArray(user->eta,&eta_p);
    283:   VecGetArray(user->logcv,&logcv_p);
    284:   VecGetArray(user->logcv2,&logcv2_p);
    285:   VecGetArray(user->DPsiv,&DPsiv_p);
    286:   VecGetArray(user->DPsieta,&DPsieta_p);

    288:   Llog(user->cv,user->logcv);

    290:   VecCopy(user->cv,user->work1);
    291:   VecScale(user->work1,-1.0);
    292:   VecShift(user->work1,1.0);
    293:   Llog(user->work1,user->logcv2);

    295:   for (i=0; i<n; i++)
    296:   {
    297:     DPsiv_p[i] = (eta_p[i]-1.0)*(eta_p[i]-1.0)*(eta_p[i]+1.0)*(eta_p[i]+1.0)*(Evf + logcv_p[i] - logcv2_p[i]) - 2.0*A*(cv_p[i] - cv0)*eta_p[i]*(eta_p[i]+2.0)*(eta_p[i]-1.0)*(eta_p[i]-1.0) + 2.0*B*(cv_p[i] - 1.0)*eta_p[i]*eta_p[i];

    299:     DPsieta_p[i] = 4.0*eta_p[i]*(eta_p[i]-1.0)*(eta_p[i]+1.0)*(Evf*cv_p[i] + cv_p[i]*logcv_p[i] + (1.0-cv_p[i])*logcv2_p[i]) - A*(cv_p[i] - cv0)*(cv_p[i] - cv0)*(4.0*eta_p[i]*eta_p[i]*eta_p[i] - 6.0*eta_p[i] + 2.0) + 2.0*B*(cv_p[i]-1.0)*(cv_p[i]-1.0)*eta_p[i];

    301:   }

    303:   VecRestoreArray(user->cv,&cv_p);
    304:   VecRestoreArray(user->eta,&eta_p);
    305:   VecGetArray(user->logcv,&logcv_p);
    306:   VecGetArray(user->logcv2,&logcv2_p);
    307:   VecRestoreArray(user->DPsiv,&DPsiv_p);
    308:   VecRestoreArray(user->DPsieta,&DPsieta_p);
    309:   return(0);

    311: }

    316: PetscErrorCode Llog(Vec X, Vec Y)
    317: {
    319:   PetscScalar    *x,*y;
    320:   PetscInt       n,i;

    323:   VecGetArray(X,&x);
    324:   VecGetArray(Y,&y);
    325:   VecGetLocalSize(X,&n);
    326:   for (i=0; i<n; i++) {
    327:     if (x[i] < 1.0e-12) y[i] = log(1.0e-12);
    328:     else y[i] = log(x[i]);
    329:   }
    330:   return(0);
    331: }

    335: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
    336: {

    339:   PetscInt          n,i,Mda;
    340:   PetscScalar       *xx,*cv_p,*wv_p,*eta_p;
    341:   PetscViewer       view_out;

    343:   /* needed for the void growth case */
    344:   PetscScalar       xmid,cv_v=1.0,cv_m=user->Sv*user->cv0,eta_v=1.0,eta_m=0.0,h,lambda;
    345:   PetscInt          nele,nen,idx[2];
    346:   const PetscInt    *ele;
    347:   PetscScalar       x[2];
    348:   Vec               coords;
    349:   const PetscScalar *_coords;
    350:   PetscScalar       xwidth = user->xmax - user->xmin;

    353:   VecGetLocalSize(X,&n);

    355:   DMDAGetInfo(user->da2,NULL,&Mda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
    356:   DMGetCoordinatesLocal(user->da2,&coords);
    357:   VecGetArrayRead(coords,&_coords);

    359:   if (user->periodic) h = (user->xmax-user->xmin)/Mda;
    360:   else                h = (user->xmax-user->xmin)/(Mda-1.0);

    362:   xmid   = (user->xmax + user->xmin)/2.0;
    363:   lambda = 4.0*h;

    365:   DMDAGetElements(user->da2,&nele,&nen,&ele);
    366:   for (i=0; i < nele; i++) {
    367:     idx[0] = ele[2*i]; idx[1] = ele[2*i+1];

    369:     x[0] = _coords[idx[0]];
    370:     x[1] = _coords[idx[1]];

    373:     PetscInt    k;
    374:     PetscScalar vals_DDcv[2],vals_cv[2],vals_eta[2],s,hhr,r;
    375:     for (k=0; k < 2; k++) {
    376:       s = PetscAbs(x[k] - xmid);
    377:       if (s <= xwidth*(5.0/64.0)) {
    378:         vals_cv[k]   = cv_v;
    379:         vals_eta[k]  = eta_v;
    380:         vals_DDcv[k] = 0.0;
    381:       } else if (s> xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0)) {
    382:         /*r = (s - xwidth*(6.0/64.0))/(0.5*lambda);*/
    383:         r            = (s - xwidth*(6.0/64.0))/(xwidth/64.0);
    384:         hhr          = 0.25*(-r*r*r + 3*r + 2);
    385:         vals_cv[k]   = cv_m + (1.0 - hhr)*(cv_v - cv_m);
    386:         vals_eta[k]  = eta_m + (1.0 - hhr)*(eta_v - eta_m);
    387:         vals_DDcv[k] = (cv_v - cv_m)*r*6.0/(lambda*lambda);
    388:       } else {
    389:         vals_cv[k]   = cv_m;
    390:         vals_eta[k]  = eta_m;
    391:         vals_DDcv[k] = 0.0;
    392:       }
    393:     }

    395:     VecSetValuesLocal(user->cv,2,idx,vals_cv,INSERT_VALUES);
    396:     VecSetValuesLocal(user->eta,2,idx,vals_eta,INSERT_VALUES);
    397:     VecSetValuesLocal(user->work2,2,idx,vals_DDcv,INSERT_VALUES);

    399:   }
    400:   DMDARestoreElements(user->da2,&nele,&nen,&ele);
    401:   VecRestoreArrayRead(coords,&_coords);

    403:   VecAssemblyBegin(user->cv);
    404:   VecAssemblyEnd(user->cv);
    405:   VecAssemblyBegin(user->eta);
    406:   VecAssemblyEnd(user->eta);
    407:   VecAssemblyBegin(user->work2);
    408:   VecAssemblyEnd(user->work2);

    410:   DPsi(user);
    411:   VecCopy(user->DPsiv,user->wv);
    412:   VecAXPY(user->wv,-2.0*user->kav,user->work2);

    414:   VecGetArray(X,&xx);
    415:   VecGetArray(user->wv,&wv_p);
    416:   VecGetArray(user->cv,&cv_p);
    417:   VecGetArray(user->eta,&eta_p);

    419:   for (i=0; i<n/3; i++) {
    420:     xx[3*i]  =wv_p[i];
    421:     xx[3*i+1]=cv_p[i];
    422:     xx[3*i+2]=eta_p[i];
    423:   }

    425:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
    426:   VecView(user->wv,view_out);
    427:   VecView(user->cv,view_out);
    428:   VecView(user->eta,view_out);
    429:   PetscViewerDestroy(&view_out);

    431:   VecRestoreArray(X,&xx);
    432:   VecRestoreArray(user->wv,&wv_p);
    433:   VecRestoreArray(user->cv,&cv_p);
    434:   VecRestoreArray(user->eta,&eta_p);
    435:   return(0);
    436: }

    440: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
    441: {
    443:   AppCtx         *user=(AppCtx*)ctx;

    446:   MatMultAdd(user->M,X,user->q,F);
    447:   return(0);
    448: }

    452: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
    453: {
    455:   AppCtx         *user=(AppCtx*)ctx;

    458:   *flg = SAME_NONZERO_PATTERN;
    459:   MatCopy(user->M,*J,*flg);
    460:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    461:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
    462:   return(0);
    463: }
    466: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
    467: {
    469:   PetscScalar    **l,**u;
    470:   PetscInt       xs,xm;
    471:   PetscInt       i;

    474:   DMDAVecGetArrayDOF(da,xl,&l);
    475:   DMDAVecGetArrayDOF(da,xu,&u);

    477:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

    480:   for (i=xs; i < xs+xm; i++) {
    481:     l[i][0] = -SNES_VI_INF;
    482:     l[i][1] = 0.0;
    483:     l[i][2] = 0.0;
    484:     u[i][0] = SNES_VI_INF;
    485:     u[i][1] = 1.0;
    486:     u[i][2] = 1.0;
    487:   }

    489:   DMDAVecRestoreArrayDOF(da,xl,&l);
    490:   DMDAVecRestoreArrayDOF(da,xu,&u);
    491:   return(0);
    492: }

    496: PetscErrorCode GetParams(AppCtx *user)
    497: {
    499:   PetscBool      flg;

    502:   /* Set default parameters */
    503:   user->xmin       = 0.0; user->xmax = 128.0;
    504:   user->Mv         = 1.0;
    505:   user->L          = 1.0;
    506:   user->kaeta      = 1.0;
    507:   user->kav        = 0.5;
    508:   user->Evf        = 9.09;
    509:   user->A          = 9.09;
    510:   user->B          = 9.09;
    511:   user->cv0        = 1.13e-4;
    512:   user->Sv         = 500.0;
    513:   user->dt         = 1.0e-5;
    514:   user->T          = 1.0e-2;
    515:   user->graphics   = PETSC_TRUE;
    516:   user->periodic   = PETSC_FALSE;
    517:   user->lumpedmass = PETSC_FALSE;

    519:   PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
    520:   PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
    521:   PetscOptionsGetReal(NULL,"-T",&user->T,&flg);
    522:   PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);
    523:   PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
    524:   PetscOptionsBool("-periodic","Use periodic boundary conditions\n","None",user->periodic,&user->periodic,&flg);
    525:   PetscOptionsBool("-lumpedmass","Use lumped mass matrix\n","None",user->lumpedmass,&user->lumpedmass,&flg);
    526:   return(0);
    527: }

    532: PetscErrorCode SetUpMatrices(AppCtx *user)
    533: {
    535:   PetscInt       nele,nen,i,n;
    536:   const PetscInt *ele;
    537:   PetscScalar    dt=user->dt,h;

    539:   PetscInt    idx[2];
    540:   PetscScalar eM_0[2][2],eM_2[2][2];
    541:   Mat         M  =user->M;
    542:   Mat         M_0=user->M_0;
    543:   PetscInt    Mda;

    547:   MatGetLocalSize(M,&n,NULL);
    548:   DMDAGetInfo(user->da1,NULL,&Mda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

    550:   if (user->periodic) h = (user->xmax-user->xmin)/Mda;
    551:   else h = (user->xmax-user->xmin)/(Mda-1.0);

    553:   if (user->lumpedmass) {
    554:     eM_0[0][0] = h/2.0;
    555:     eM_0[1][1] = h/2.0;
    556:     eM_0[0][1] = eM_0[1][0] = 0.0;
    557:   } else {
    558:     eM_0[0][0]=eM_0[1][1]=h/3.0;
    559:     eM_0[0][1]=eM_0[1][0]=h/6.0;
    560:   }
    561:   eM_2[0][0]=eM_2[1][1]=1.0/h;
    562:   eM_2[0][1]=eM_2[1][0]=-1.0/h;

    564:   /* Get local element info */
    565:   DMDAGetElements(user->da1,&nele,&nen,&ele);
    566:   for (i=0; i < nele; i++) {

    568:     idx[0] = ele[2*i]; idx[1] = ele[2*i+1];

    570:     PetscInt    row,cols[4],r,row_M_0,cols2[2];
    571:     PetscScalar vals[4],vals_M_0[2],vals2[2];

    573:     for (r=0; r<2; r++) {
    574:       row_M_0    = idx[r];
    575:       vals_M_0[0]=eM_0[r][0];
    576:       vals_M_0[1]=eM_0[r][1];

    578:       MatSetValuesLocal(M_0,1,&row_M_0,2,idx,vals_M_0,ADD_VALUES);

    580:       row     = 3*idx[r];
    581:       cols[0] = 3*idx[0];     vals[0] = dt*eM_2[r][0]*user->Mv;
    582:       cols[1] = 3*idx[1];     vals[1] = dt*eM_2[r][1]*user->Mv;
    583:       cols[2] = 3*idx[0]+1;   vals[2] = eM_0[r][0];
    584:       cols[3] = 3*idx[1]+1;   vals[3] = eM_0[r][1];

    586:       /* Insert values in matrix M for 1st dof */
    587:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);

    589:       row     = 3*idx[r]+1;
    590:       cols[0] = 3*idx[0];     vals[0] = -eM_0[r][0];
    591:       cols[1] = 3*idx[1];     vals[1] = -eM_0[r][1];
    592:       cols[2] = 3*idx[0]+1;   vals[2] = 2.0*user->kav*eM_2[r][0];
    593:       cols[3] = 3*idx[1]+1;   vals[3] = 2.0*user->kav*eM_2[r][1];

    595:       /* Insert values in matrix M for 2nd dof */
    596:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);

    599:       row      = 3*idx[r]+2;
    600:       cols2[0] = 3*idx[0]+2;   vals2[0] = eM_0[r][0] + user->dt*2.0*user->L*user->kaeta*eM_2[r][0];
    601:       cols2[1] = 3*idx[1]+2;   vals2[1] = eM_0[r][1] + user->dt*2.0*user->L*user->kaeta*eM_2[r][1];

    603:       MatSetValuesLocal(M,1,&row,2,cols2,vals2,ADD_VALUES);
    604:     }
    605:   }

    607:   MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
    608:   MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);

    610:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
    611:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

    613:   DMDARestoreElements(user->da1,&nele,&nen,&ele);
    614:   return(0);
    615: }

    619: PetscErrorCode CheckRedundancy(SNES snes, IS act, IS *outact, DM da)
    620: {
    622:   PetscScalar    **uin,**uout;
    623:   Vec            UIN, UOUT;
    624:   PetscInt       xs,xm,*outindex;
    625:   const PetscInt *index;
    626:   PetscInt       k,i,l,n,M,cnt=0;

    629:   DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
    630:   DMGetGlobalVector(da,&UIN);
    631:   VecSet(UIN,0.0);
    632:   DMGetLocalVector(da,&UOUT);
    633:   DMDAVecGetArrayDOF(da,UIN,&uin);
    634:   ISGetIndices(act,&index);
    635:   ISGetLocalSize(act,&n);
    636:   for (k=0; k<n; k++) {
    637:     l = index[k]%5;
    638:     i = index[k]/5;

    640:     uin[i][l]=1.0;
    641:   }
    642:   printf("Number of active constraints before applying redundancy %d\n",n);
    643:   ISRestoreIndices(act,&index);
    644:   DMDAVecRestoreArrayDOF(da,UIN,&uin);
    645:   DMGlobalToLocalBegin(da,UIN,INSERT_VALUES,UOUT);
    646:   DMGlobalToLocalEnd(da,UIN,INSERT_VALUES,UOUT);
    647:   DMDAVecGetArrayDOF(da,UOUT,&uout);

    649:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

    651:   for (i=xs; i < xs+xm;i++) {
    652:     if (uout[i-1][1] && uout[i][1] && uout[i+1][1]) uout[i][0] = 1.0;
    653:     if (uout[i-1][3] && uout[i][3] && uout[i+1][3]) uout[i][2] = 1.0;
    654:   }

    656:   for (i=xs; i < xs+xm; i++) {
    657:     for (l=0; l < 5; l++) {
    658:       if (uout[i][l]) cnt++;
    659:     }
    660:   }

    662:   printf("Number of active constraints after applying redundancy %d\n",cnt);

    664:   PetscMalloc(cnt*sizeof(PetscInt),&outindex);
    665:   cnt  = 0;

    667:   for (i=xs; i < xs+xm;i++) {
    668:     for (l=0;l<5;l++) {
    669:       if (uout[i][l]) outindex[cnt++] = 5*(i)+l;
    670:     }
    671:   }

    674:   ISCreateGeneral(PETSC_COMM_WORLD,cnt,outindex,PETSC_OWN_POINTER,outact);
    675:   DMDAVecRestoreArrayDOF(da,UOUT,&uout);
    676:   DMRestoreGlobalVector(da,&UIN);
    677:   DMRestoreLocalVector(da,&UOUT);
    678:   return(0);
    679: }