Actual source code: ex30.cxx

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: static const char help[] = "Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods and uses MOAB.\n";
  2: /*
  3:    u_t - alpha u_xx = A + u^2 v - (B+1) u
  4:    v_t - alpha v_xx = B u - u^2 v
  5:    0 < x < 1;
  6:    A = 1, B = 3, alpha = 1/50

  8:    Initial conditions:
  9:    u(x,0) = 1 + sin(2 pi x)
 10:    v(x,0) = 3

 12:    Boundary conditions:
 13:    u(0,t) = u(1,t) = 1
 14:    v(0,t) = v(1,t) = 3
 15: */

 17: // PETSc includes:
 18: #include <petscts.h>
 19: #include <petscdmmoab.h>

 21: // MOAB includes:
 22: #if defined (PETSC_HAVE_MOAB)
 23: #  include <moab/Core.hpp>
 24: #  include <moab/ReadUtilIface.hpp>
 25: #  include <MBTagConventions.hpp>
 26: #else
 27: #error You must have MOAB for this example. Reconfigure using --download-moab
 28: #endif

 30: typedef struct {
 31:   PetscScalar u,v;
 32: } Field;

 34: typedef struct _User *User;
 35: struct _User {
 36:   PetscReal A,B;                /* Reaction coefficients */
 37:   PetscReal alpha;              /* Diffusion coefficient */
 38:   PetscReal uleft,uright;       /* Dirichlet boundary conditions */
 39:   PetscReal vleft,vright;       /* Dirichlet boundary conditions */
 40:   PetscInt  npts;               /* Number of mesh points */
 41:   PetscBool io;

 43:   moab::ParallelComm *pcomm;
 44:   moab::Interface *mbint;
 45:   moab::Range *owned_vertexes;
 46:   moab::Range *owned_edges;
 47:   moab::Range *all_vertexes;
 48:   moab::Range *shared_vertexes;
 49:   moab::Tag    unknowns_tag;
 50:   PetscInt     unknowns_tag_size;
 51:   moab::Tag    id_tag;
 52: };

 54: static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
 55: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
 56: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);

 58: PetscErrorCode create_app_data(_User& user);
 59: PetscErrorCode destroy_app_data(_User& user);
 60: PetscErrorCode create_matrix(_User &user, DM &dm, Mat *J);

 62: /****************
 63:  *              *
 64:  *     MAIN     *
 65:  *              *
 66:  ****************/
 69: int main(int argc,char **argv)
 70: {
 71:   TS                ts;         /* nonlinear solver */
 72:   Vec               X;          /* solution, residual vectors */
 73:   Mat               J;          /* Jacobian matrix */
 74:   PetscInt          steps,maxsteps,mx;
 75:   PetscErrorCode    ierr;
 76:   PetscReal         ftime,hx,dt;
 77:   _User             user;       /* user-defined work context */
 78:   TSConvergedReason reason;
 79:   DM                dm;
 80:   moab::ErrorCode   merr;

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

 84:   /* Fill in the user defined work context (creates mesh, solution field on mesh) */
 85:   create_app_data(user);

 87:   /* Create the DM to manage the mesh */
 88:   DMMoabCreateMoab(user.pcomm->comm(),user.mbint,user.pcomm,&user.id_tag,user.owned_vertexes,&dm);
 89:   DMSetUp(dm);

 91:   /* Create the solution vector: */
 92:   DMMoabCreateVector(dm,user.unknowns_tag,user.owned_vertexes,PETSC_TRUE,PETSC_FALSE,&X);
 93:   /* Create the matrix */
 94:   create_matrix(user, dm, &J);

 96:   /* Create timestepping solver context */
 97:   TSCreate(PETSC_COMM_WORLD,&ts);
 98:   TSSetType(ts,TSARKIMEX);

100:   TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);
101:   TSSetIFunction(ts,NULL,FormIFunction,&user);
102:   TSSetIJacobian(ts,J,J,FormIJacobian,&user);
103:   TSSetDM(ts,dm);

105:   ftime = 10.0;
106:   maxsteps = 10000;
107:   TSSetDuration(ts,maxsteps,ftime);

109:   /* Set initial conditions */
110:   TSSetSolution(ts,X);
111:   VecGetSize(X,&mx);
112:   hx = 1.0/(PetscReal)(mx/2-1);
113:   dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */
114:   TSSetInitialTimeStep(ts,0.0,dt);

116:   /* Set runtime options */
117:   TSSetFromOptions(ts);

119:   /* Solve nonlinear system */
120:   TSSolve(ts,X);
121:   TSGetSolveTime(ts,&ftime);
122:   TSGetTimeStepNumber(ts,&steps);
123:   TSGetConvergedReason(ts,&reason);
124:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);

126:   /* Write out the final mesh */
127:   if (user.io) {
128: #ifdef MOAB_HDF5_H
129:     merr = user.mbint->write_file("ex30.h5m");MBERRNM(merr);
130: #else
131:     merr = user.mbint->write_file("ex30.vtk");MBERRNM(merr);
132: #endif
133:   }

135:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);

137:   /* Free work space. */
138:   MatDestroy(&J);
139:   VecDestroy(&X);
140:   TSDestroy(&ts);
141:   DMDestroy(&dm);

143:   /* Free all MOAB related resources */
144:   destroy_app_data(user);

146:   PetscFinalize();
147:   return 0;
148: }

152: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
153: {
154:   User            user = (User)ptr;
155:   PetscInt        i;
156:   Field           *x,*xdot,*f;
157:   PetscReal       hx;
158:   PetscErrorCode  ierr;
159:   moab::ErrorCode merr;
160:   int             *vertex_ids;

163:   hx = 1.0/(PetscReal)(user->npts-1);

165:   // Get connectivity information:
166:   int verts_per_entity;
167:   int count,num_edges;
168:   moab::EntityHandle *connect;
169:   moab::Range::iterator iter = user->owned_edges->begin();
170:   merr = user->mbint->connect_iterate(iter,user->owned_edges->end(),connect,verts_per_entity,num_edges);MBERRNM(merr);

172:   // Get tag data:
173:   moab::Tag x_tag;
174:   DMMoabGetVecTag(X,&x_tag);
175:   merr = user->mbint->tag_iterate(x_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
176:                                   count,reinterpret_cast<void*&>(x));MBERRNM(merr);

178:   moab::Tag xdot_tag;
179:   DMMoabGetVecTag(Xdot,&xdot_tag);
180:   merr = user->mbint->tag_iterate(xdot_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
181:                                   count,reinterpret_cast<void*&>(xdot));MBERRNM(merr);

183:   moab::Tag f_tag;
184:   DMMoabGetVecTag(F,&f_tag);
185:   merr = user->mbint->tag_iterate(f_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
186:                                   count,reinterpret_cast<void*&>(f));MBERRNM(merr);

188:   // Get the global IDs on all vertexes:
189:   moab::Tag id_tag;
190:   user->mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);
191:   merr = user->mbint->tag_iterate(id_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
192:                                     count,reinterpret_cast<void*&>(vertex_ids));MBERRNM(merr);

194:   // Exchange tags that are needed for assembly:
195:   merr = user->pcomm->exchange_tags(x_tag,*user->all_vertexes);MBERRNM(merr);
196:   merr = user->pcomm->exchange_tags(xdot_tag,*user->all_vertexes);MBERRNM(merr);

198:   // Compute f:
199:   const Field zero_field = {0.0,0.0};
200:   std::fill(f,f+user->all_vertexes->size(),zero_field);

202:   const moab::EntityHandle first_vertex = *user->all_vertexes->begin();

204:   for (i = 0; i < num_edges; i++) {
205:     const moab::EntityHandle idx_left  = connect[2*i]-first_vertex;
206:     const moab::EntityHandle idx_right = connect[2*i+1]-first_vertex;

208:     const int id_left  = vertex_ids[idx_left ];
209:     const int id_right = vertex_ids[idx_right];

211:     if (id_left == 0) {
212:       // Apply left BC
213:       f[idx_left].u += hx * (x[idx_left].u - user->uleft);
214:       f[idx_left].v += hx * (x[idx_left].v - user->vleft);
215:     } else {
216:       f[idx_left].u += hx * xdot[idx_left].u - user->alpha*(-2*x[idx_left].u + x[idx_right].u)/hx;
217:       f[idx_left].v += hx * xdot[idx_left].v - user->alpha*(-2*x[idx_left].v + x[idx_right].v)/hx;
218:     }

220:     if (id_right == user->npts-1) {
221:       // Apply right BC
222:       f[idx_right].u += hx * (x[idx_right].u - user->uright);
223:       f[idx_right].v += hx * (x[idx_right].v - user->vright);
224:     } else {
225:       f[idx_right].u -= user->alpha*x[idx_left].u/hx;
226:       f[idx_right].v -= user->alpha*x[idx_left].v/hx;
227:     }
228:   }

230:   // Add tags on shared vertexes:
231:   merr = user->pcomm->reduce_tags(f_tag,MPI_SUM,*user->shared_vertexes);MBERRNM(merr);
232:   return(0);
233: }

237: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
238: {
239:   User            user = (User)ptr;
240:   PetscReal       hx;
241:   Field           *x,*f;
242:   PetscErrorCode  ierr;
243:   moab::Tag       id_tag;
244:   moab::ErrorCode merr;

247:   hx = 1.0/(PetscReal)(user->npts-1);

249:   merr = user->mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);MBERRNM(merr);

251:   /* Get pointers to vector data */
252:   VecGetArray(X,reinterpret_cast<PetscScalar**>(&x));
253:   VecGetArray(F,reinterpret_cast<PetscScalar**>(&f));

255:   /* Compute function over the locally owned part of the grid */
256:   const moab::EntityHandle first_vertex = *user->owned_vertexes->begin();
257:   moab::Range::iterator iter;
258:   for(iter = user->owned_vertexes->begin(); iter != user->owned_vertexes->end(); iter++) {
259:     moab::EntityHandle i = *iter - first_vertex;
260:     PetscScalar u = x[i].u,v = x[i].v;
261:     f[i].u = hx*(user->A + u*u*v - (user->B+1)*u);
262:     f[i].v = hx*(user->B*u - u*u*v);
263:   }

265:   /* Restore vectors */
266:   VecRestoreArray(X,reinterpret_cast<PetscScalar**>(&x));
267:   VecRestoreArray(F,reinterpret_cast<PetscScalar**>(&f));
268:   return(0);
269: }

271: /* --------------------------------------------------------------------- */
272: /*
273:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
274: */
277: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ptr)
278: {
279:   User            user = (User)ptr;
280:   PetscErrorCode  ierr;
281:   PetscInt        i;
282:   PetscReal       hx;
283:   moab::Tag       id_tag;
284:   moab::ErrorCode merr;

287:   hx = 1.0/(PetscReal)(user->npts-1);

289:   merr = user->mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);MBERRNM(merr);

291:   /* Compute function over the locally owned part of the grid */
292:   moab::Range::iterator iter;
293:   for(iter = user->owned_vertexes->begin(); iter != user->owned_vertexes->end(); iter++) {
294:     merr = user->mbint->tag_get_data(id_tag,&(*iter),1,&i);MBERRNM(merr);

296:     if (i == 0 || i == user->npts-1) {
297:       // Boundary conditions...
298:       const PetscInt row = i,col = i;
299:       const PetscScalar vals[2][2] = {{hx,0},{0,hx}};
300:       MatSetValuesBlocked(Jpre,1,&row,1,&col,&vals[0][0],INSERT_VALUES);
301:     } else {
302:       //
303:       const PetscInt row = i,col[] = {i-1,i,i+1};
304:       const PetscScalar dxxL = -user->alpha/hx,dxx0 = 2.*user->alpha/hx,dxxR = -user->alpha/hx;
305:       const PetscScalar vals[2][3][2] = {{{dxxL,0},{a*hx+dxx0,0},{dxxR,0}},
306:                                          {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}};
307:       MatSetValuesBlocked(Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);
308:     }
309:   }

311:   MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
312:   MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
313:   if (J != Jpre) {
314:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
315:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
316:   }
317:   return(0);
318: }


321: /********************************
322:  *                              *
323:  *     initialize_moab_mesh     *
324:  *                              *
325:  ********************************/

329: PetscErrorCode initialize_moab_mesh(moab::ParallelComm* pcomm,int npts,int nghost,moab::Tag &unknowns_tag,PetscInt &unknowns_tag_size,moab::Tag &id_tag)
330: {
331:   moab::ErrorCode merr;
332:   PetscInt num_procs;
333:   PetscInt rank;

336:   MPI_Comm_size( pcomm->comm(),&num_procs );
337:   MPI_Comm_rank( pcomm->comm(),&rank );

339:   // Begin with some error checking:
340:   if(npts < 2) {
341:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"npts must be >= 2");
342:   }

344:   if(num_procs >= npts) {
345:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Num procs must be < npts");
346:   }

348:   // No errors,proceed with building the mesh:
349:   moab::Interface *mbint = pcomm->get_moab();
350:   moab::ReadUtilIface* readMeshIface;
351:   mbint->query_interface(readMeshIface);

353:   // Determine which elements (cells) this process owns:
354:   const PetscInt nele = npts-1;
355:   PetscInt my_nele = nele / num_procs; // Number of elements owned by this process
356:   const PetscInt extra = nele % num_procs;
357:   PetscInt vstart = rank * my_nele;    // The starting element for this process
358:   if(rank < extra) my_nele++;
359:   vstart += std::min(extra, rank);
360:   const PetscInt my_npts = my_nele + 1;

362:   // Create the local portion of the mesh:
363:   // Create vertexes and set the coodinate of each vertex:
364:   moab::EntityHandle vertex_first;
365:   std::vector<double*> vertex_coords;
366:   const int sequence_size = (my_nele + 2) + 1;
367:   merr = readMeshIface->get_node_coords(3,my_npts,1,vertex_first,vertex_coords,sequence_size);MBERRNM(merr);

369:   const double xs = 0.0, xe = 1.0;
370:   const double dx = (xe - xs) / nele;
371:   for (int i = 0; i < my_npts; i++) {
372:     vertex_coords[0][i] = (i+vstart)*dx;
373:     vertex_coords[1][i] = vertex_coords[2][i] = 0.0;
374:   }

376:   moab::Range owned_vertexes;
377:   merr = mbint->get_entities_by_type(0,moab::MBVERTEX,owned_vertexes);MBERRNM(merr);

379:   // Create elements between mesh points. This is done so that VisIt
380:   // will interpret the output as a mesh that can be plotted...
381:   moab::EntityHandle edge_first;
382:   moab::EntityHandle *connectivity = 0;
383:   merr = readMeshIface->get_element_connect(my_nele,2,moab::MBEDGE,1,edge_first,connectivity);MBERRNM(merr);

385:   for (int i = 0; i < my_nele; i+=1) {
386:     connectivity[2*i]   = vertex_first + i;
387:     connectivity[2*i+1] = vertex_first + (i+1);
388:   }

390:   merr = readMeshIface->update_adjacencies(edge_first,my_nele,2,connectivity);MBERRNM(merr);

392:   // Set tags on all of the vertexes...
393:   // Create tag handle to represent the unknowns, u and v and
394:   // initialize them. We will create a single tag whose type is an
395:   // array of two doubles and whose name is "unknowns"
396:   Field default_value = {0.0,0.0};
397:   unknowns_tag_size = sizeof(Field)/sizeof(PetscScalar);
398:   merr = mbint->tag_get_handle("unknowns",unknowns_tag_size,moab::MB_TYPE_DOUBLE,unknowns_tag,
399:                               moab::MB_TAG_DENSE | moab::MB_TAG_CREAT,&default_value);MBERRNM(merr);

401:   // Apply the "unknowns" tag to the vertexes with some initial value...
402:   std::vector<Field> tag_data(my_npts);
403:   for (int i = 0; i < my_npts; i++) {
404:     tag_data[i].u = 1+sin(2*PETSC_PI*vertex_coords[0][i]);
405:     tag_data[i].v = 3;
406:   }

408:   merr = mbint->tag_set_data(unknowns_tag,owned_vertexes,tag_data.data());MBERRNM(merr);

410:   // Get the global ID tag. The global ID tag is applied to each
411:   // vertex. It acts as an global identifier which MOAB uses to
412:   // assemble the individual pieces of the mesh:
413:   merr = mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);MBERRNM(merr);

415:   std::vector<int> global_ids(my_npts);
416:   for (int i = 0; i < my_npts; i++) {
417:     global_ids[i] = i+vstart;
418:   }

420:   merr = mbint->tag_set_data(id_tag,owned_vertexes,global_ids.data());MBERRNM(merr);

422:   moab::Range owned_edges;
423:   merr = mbint->get_entities_by_type(0,moab::MBEDGE,owned_edges);MBERRNM(merr);
424:   merr = pcomm->resolve_shared_ents(0,owned_edges,1,0);MBERRNM(merr);

426:   // Reassign global IDs on all entities.
427:   merr = pcomm->assign_global_ids(0,1,0,false);MBERRNM(merr);
428:   merr = pcomm->exchange_ghost_cells( 1,0,nghost,0,true);MBERRNM(merr);

430:   // Everything is set up, now just do a tag exchange to update tags
431:   // on all of the shared vertexes:
432:   merr = pcomm->exchange_tags(id_tag,owned_vertexes);MBERRNM(merr);
433:   return(0);
434: }


439: PetscErrorCode create_app_data(_User& user)
440: {
442:   moab::ErrorCode merr;


446:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Advection-reaction options","ex30.cxx"); {
447:     user.A      = 1;
448:     user.B      = 3;
449:     user.alpha  = 0.02;
450:     user.uleft  = 1;
451:     user.uright = 1;
452:     user.vleft  = 3;
453:     user.vright = 3;
454:     user.npts   = 11;
455:     user.io     = PETSC_FALSE;
456:     PetscOptionsReal("-A","Reaction rate","ex30.cxx",user.A,&user.A,NULL);
457:     PetscOptionsReal("-B","Reaction rate","ex30.cxx",user.B,&user.B,NULL);
458:     PetscOptionsReal("-alpha","Diffusion coefficient","ex30.cxx",user.alpha,&user.alpha,NULL);
459:     PetscOptionsReal("-uleft","Dirichlet boundary condition","ex30.cxx",user.uleft,&user.uleft,NULL);
460:     PetscOptionsReal("-uright","Dirichlet boundary condition","ex30.cxx",user.uright,&user.uright,NULL);
461:     PetscOptionsReal("-vleft","Dirichlet boundary condition","ex30.cxx",user.vleft,&user.vleft,NULL);
462:     PetscOptionsReal("-vright","Dirichlet boundary condition","ex30.cxx",user.vright,&user.vright,NULL);
463:     PetscOptionsInt("-npts","Number of mesh points","ex30.cxx",user.npts,&user.npts,NULL);
464:     PetscOptionsBool("-io", "Write out the solution and mesh data", "ex30.cxx", user.io, &user.io, NULL);
465:   } PetscOptionsEnd();

467:   user.mbint = new moab::Core;
468:   user.pcomm = new moab::ParallelComm(user.mbint,PETSC_COMM_WORLD);

470:   initialize_moab_mesh(user.pcomm,user.npts,1,user.unknowns_tag,user.unknowns_tag_size,user.id_tag);  // creates moab mesh and field of unknowns on vertices

472:   // Set the edge/vertex range once so we don't have to do it
473:   // again. To do this we get all of the edges/vertexes then filter so
474:   // we have only the owned entities in the ranges:
475:   user.owned_edges = new moab::Range;
476:   merr = user.mbint->get_entities_by_type(0,moab::MBEDGE,*user.owned_edges);MBERRNM(merr);

478:   user.owned_vertexes = new moab::Range;
479:   merr = user.mbint->get_entities_by_type(0,moab::MBVERTEX,*user.owned_vertexes);MBERRNM(merr);
480:   user.all_vertexes = new moab::Range(*user.owned_vertexes);

482:   user.shared_vertexes = new moab::Range;
483:   merr = user.mbint->get_entities_by_type(0,moab::MBVERTEX,*user.shared_vertexes);MBERRNM(merr);

485:   merr = user.pcomm->filter_pstatus(*user.owned_edges,PSTATUS_NOT_OWNED,PSTATUS_NOT);MBERRNM(merr);

487:   merr = user.pcomm->filter_pstatus(*user.owned_vertexes,PSTATUS_NOT_OWNED,PSTATUS_NOT);MBERRNM(merr);

489:   merr = user.pcomm->filter_pstatus(*user.shared_vertexes,PSTATUS_SHARED,PSTATUS_OR);MBERRNM(merr);

491:   // Do some error checking...make sure that tag_data is in a single
492:   // sequence:
493:   int count;
494:   void *data;
495:   moab::Tag tag;
496:   merr = user.mbint->tag_get_handle("unknowns",tag);MBERRNM(merr);
497:   merr = user.mbint->tag_iterate(tag,user.all_vertexes->begin(),user.all_vertexes->end(),
498:                                  count,data);MBERRNM(merr);
499:   if((unsigned int) count != user.all_vertexes->size()) {
500:     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Tag data not laid out contiguously %i %i",
501:              count,user.all_vertexes->size());
502:   }
503:   return(0);
504: }


509: PetscErrorCode create_matrix(_User &user, DM &dm, Mat *J)
510: {
512:   moab::ErrorCode merr;
513:   moab::Tag ltog_tag;
514:   const moab::Range *range;

517:   DMMoabGetLocalToGlobalTag(dm,&ltog_tag);
518:   DMMoabGetLocalVertices(dm,&range,NULL);

520:   // Create the matrix:
521:   MatCreateBAIJ(user.pcomm->comm(),2,2*range->size(),2*range->size(),
522:                          PETSC_DECIDE,PETSC_DECIDE,3, NULL, 3, NULL, J);

524:   // Set local to global numbering using the ltog_tag:
525:   ISLocalToGlobalMapping ltog;
526:   PetscInt               *gindices = new PetscInt[range->size()];
527:   merr = user.pcomm->get_moab()->tag_get_data(ltog_tag, *range, gindices);MBERRNM(merr);

529:   ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 2, range->size(), gindices,PETSC_COPY_VALUES, &ltog);
530:   MatSetLocalToGlobalMapping(*J,ltog,ltog);

532:   // Clean up:
533:   ISLocalToGlobalMappingDestroy(&ltog);
534:   delete [] gindices;
535:   return(0);
536: }

540: PetscErrorCode destroy_app_data(_User& user)
541: {

544:   delete user.owned_edges;
545:   delete user.owned_vertexes;
546:   delete user.all_vertexes;
547:   delete user.shared_vertexes;
548:   delete user.pcomm;
549:   delete user.mbint;

551:   return(0);
552: }