Actual source code: ex30.cxx

petsc-3.4.4 2014-03-13
  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 */

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

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

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

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

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

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

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

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

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

100:   TSSetRHSFunction(ts,PETSC_NULL,FormRHSFunction,&user);
101:   TSSetIFunction(ts,PETSC_NULL,FormIFunction,&user);
102:   TSSetIJacobian(ts,J,J,FormIJacobian,&user);

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

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

117:   // /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:   //    Set runtime options
119:   //  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120:   TSSetFromOptions(ts);

122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123:      Solve nonlinear system
124:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125:   TSSolve(ts,X);
126:   TSGetSolveTime(ts,&ftime);
127:   TSGetTimeStepNumber(ts,&steps);
128:   TSGetConvergedReason(ts,&reason);
129:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %G after %D steps\n",TSConvergedReasons[reason],ftime,steps);


132:   // Print the solution:
133:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:      Write out the final mesh
137:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138:   merr = user.mbint->write_file("ex30-final.h5m");MBERRNM(merr);

140:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141:      Free work space.
142:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

144:   // Free all PETSc related resources:
145:   MatDestroy(&J);
146:   VecDestroy(&X);
147:   TSDestroy(&ts);
148:   DMDestroy(&dm);

150:   // Free all MOAB related resources:
151:   destroy_app_data(user);

153:   PetscFinalize();
154:   return 0;
155: }

159: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
160: {
161:   User            user = (User)ptr;
162:   PetscInt        i;
163:   Field           *x,*xdot,*f;
164:   PetscReal       hx;
165:   PetscErrorCode  ierr;
166:   moab::ErrorCode merr;
167:   int             *vertex_ids;

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

172:   // Get connectivity information:
173:   int verts_per_entity;
174:   int count,num_edges;
175:   moab::EntityHandle *connect;
176:   moab::Range::iterator iter = user->owned_edges->begin();
177:   merr = user->mbint->connect_iterate(iter,user->owned_edges->end(),connect,verts_per_entity,num_edges);MBERRNM(merr);

179:   // Get tag data:
180:   moab::Tag x_tag;
181:   DMMoabGetVecTag(X,&x_tag);
182:   merr = user->mbint->tag_iterate(x_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
183:                                   count,reinterpret_cast<void*&>(x));MBERRNM(merr);

185:   moab::Tag xdot_tag;
186:   DMMoabGetVecTag(Xdot,&xdot_tag);
187:   merr = user->mbint->tag_iterate(xdot_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
188:                                   count,reinterpret_cast<void*&>(xdot));MBERRNM(merr);

190:   moab::Tag f_tag;
191:   DMMoabGetVecTag(F,&f_tag);
192:   merr = user->mbint->tag_iterate(f_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
193:                                   count,reinterpret_cast<void*&>(f));MBERRNM(merr);

195:   // Get the global IDs on all vertexes:
196:   moab::Tag id_tag;
197:   user->mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);
198:   merr = user->mbint->tag_iterate(id_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
199:                                     count,reinterpret_cast<void*&>(vertex_ids));MBERRNM(merr);

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

205:   // Compute f:
206:   const Field zero_field = {0.0,0.0};
207:   std::fill(f,f+user->all_vertexes->size(),zero_field);

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

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

215:     const int id_left  = vertex_ids[idx_left ];
216:     const int id_right = vertex_ids[idx_right];

218:     if (id_left == 0) {
219:       // Apply left BC
220:       f[idx_left].u += hx * (x[idx_left].u - user->uleft);
221:       f[idx_left].v += hx * (x[idx_left].v - user->vleft);
222:     } else {
223:       f[idx_left].u += hx * xdot[idx_left].u - user->alpha*(-2*x[idx_left].u + x[idx_right].u)/hx;
224:       f[idx_left].v += hx * xdot[idx_left].v - user->alpha*(-2*x[idx_left].v + x[idx_right].v)/hx;
225:     }

227:     if (id_right == user->npts-1) {
228:       // Apply right BC
229:       f[idx_right].u += hx * (x[idx_right].u - user->uright);
230:       f[idx_right].v += hx * (x[idx_right].v - user->vright);
231:     } else {
232:       f[idx_right].u -= user->alpha*x[idx_left].u/hx;
233:       f[idx_right].v -= user->alpha*x[idx_left].v/hx;
234:     }
235:   }

237:   // Add tags on shared vertexes:
238:   merr = user->pcomm->reduce_tags(f_tag,MPI_SUM,*user->shared_vertexes);MBERRNM(merr);

240:   // Print vectors for debugging:
241:   // PetscPrintf(PETSC_COMM_WORLD, "X\n");
242:   // VecView(X,PETSC_VIEWER_STDOUT_WORLD);

244:   // PetscPrintf(PETSC_COMM_WORLD, "Xdot\n");
245:   // VecView(Xdot,PETSC_VIEWER_STDOUT_WORLD);

247:   // PetscPrintf(PETSC_COMM_WORLD, "F\n");
248:   // VecView(F,PETSC_VIEWER_STDOUT_WORLD);

250:   return(0);
251: }

255: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
256: {
257:   User            user = (User)ptr;
258:   PetscReal       hx;
259:   Field           *x,*f;
260:   PetscErrorCode  ierr;
261:   moab::Tag       id_tag;
262:   moab::ErrorCode merr;

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

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

269:   /* Get pointers to vector data */
270:   VecGetArray(X,reinterpret_cast<PetscScalar**>(&x));
271:   VecGetArray(F,reinterpret_cast<PetscScalar**>(&f));

273:   /* Compute function over the locally owned part of the grid */
274:   const moab::EntityHandle first_vertex = *user->owned_vertexes->begin();
275:   moab::Range::iterator iter;
276:   for(iter = user->owned_vertexes->begin(); iter != user->owned_vertexes->end(); iter++) {
277:     moab::EntityHandle i = *iter - first_vertex;
278:     PetscScalar u = x[i].u,v = x[i].v;
279:     f[i].u = hx*(user->A + u*u*v - (user->B+1)*u);
280:     f[i].v = hx*(user->B*u - u*u*v);
281:   }

283:   /* Restore vectors */
284:   VecRestoreArray(X,reinterpret_cast<PetscScalar**>(&x));
285:   VecRestoreArray(F,reinterpret_cast<PetscScalar**>(&f));

287:   // Print vectors for debugging:
288:   /* PetscPrintf(PETSC_COMM_WORLD,"RHS"); */
289:   /* VecView(F,PETSC_VIEWER_STDOUT_WORLD); */

291:   return(0);
292: }

294: /* --------------------------------------------------------------------- */
295: /*
296:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
297: */
300: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *J,Mat *Jpre,MatStructure *str,void *ptr)
301: {
302:   User            user = (User)ptr;
303:   PetscErrorCode  ierr;
304:   PetscInt        i;
305:   PetscReal       hx;
306:   moab::Tag       id_tag;
307:   moab::ErrorCode merr;

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

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

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

319:     if (i == 0 || i == user->npts-1) {
320:       // Boundary conditions...
321:       const PetscInt row = i,col = i;
322:       const PetscScalar vals[2][2] = {{hx,0},{0,hx}};
323:       MatSetValuesBlocked(*Jpre,1,&row,1,&col,&vals[0][0],INSERT_VALUES);
324:     } else {
325:       //
326:       const PetscInt row = i,col[] = {i-1,i,i+1};
327:       const PetscScalar dxxL = -user->alpha/hx,dxx0 = 2.*user->alpha/hx,dxxR = -user->alpha/hx;
328:       const PetscScalar vals[2][3][2] = {{{dxxL,0},{a*hx+dxx0,0},{dxxR,0}},
329:                                          {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}};
330:       MatSetValuesBlocked(*Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);
331:     }
332:   }

334:   MatAssemblyBegin(*Jpre,MAT_FINAL_ASSEMBLY);
335:   MatAssemblyEnd(*Jpre,MAT_FINAL_ASSEMBLY);
336:   if (*J != *Jpre) {
337:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
338:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
339:   }

341:   /* MatView(*J,PETSC_VIEWER_STDOUT_WORLD); */

343:   return(0);
344: }


347: /********************************
348:  *                              *
349:  *     initialize_moab_mesh     *
350:  *                              *
351:  ********************************/

355: PetscErrorCode initialize_moab_mesh(moab::ParallelComm* pcomm,int npts,int nghost,moab::Tag &unknowns_tag,PetscInt &unknowns_tag_size,moab::Tag &id_tag)
356: {
357:   moab::ErrorCode merr;
358:   PetscInt num_procs;
359:   PetscInt rank;

362:   MPI_Comm_size( pcomm->comm(),&num_procs );
363:   MPI_Comm_rank( pcomm->comm(),&rank );

365:   // Begin with some error checking:
366:   if(npts < 2) {
367:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"npts must be >= 2");
368:   }

370:   if(num_procs >= npts) {
371:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Num procs must be < npts");
372:   }

374:   // No errors,proceed with building the mesh:
375:   moab::Interface *mbint = pcomm->get_moab();
376:   moab::ReadUtilIface* readMeshIface;
377:   mbint->query_interface(readMeshIface);

379:   // Determine which elements (cells) this process owns:
380:   const PetscInt nele = npts-1;
381:   PetscInt my_nele = nele / num_procs; // Number of elements owned by this process
382:   const PetscInt extra = nele % num_procs;
383:   PetscInt vstart = rank * my_nele;    // The starting element for this process
384:   if(rank < extra) my_nele++;
385:   vstart += std::min(extra, rank);
386:   const PetscInt my_npts = my_nele + 1;

388:   // Create the local portion of the mesh:

390:   // Create vertexes and set the coodinate of each vertex:
391:   moab::EntityHandle vertex_first;
392:   std::vector<double*> vertex_coords;
393:   const int sequence_size = (my_nele + 2) + 1;
394:   merr = readMeshIface->get_node_coords(3,my_npts,1,vertex_first,vertex_coords,sequence_size);MBERRNM(merr);

396:   const double xs = 0.0, xe = 1.0;
397:   const double dx = (xe - xs) / nele;
398:   for (int i = 0; i < my_npts; i++) {
399:     vertex_coords[0][i] = (i+vstart)*dx;
400:     vertex_coords[1][i] = vertex_coords[2][i] = 0.0;
401:   }

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

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

412:   for (int i = 0; i < my_nele; i+=1) {
413:     connectivity[2*i]   = vertex_first + i;
414:     connectivity[2*i+1] = vertex_first + (i+1);
415:   }

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

419:   // Set tags on all of the vertexes...

421:   // Create tag handle to represent the unknowns, u and v and
422:   // initialize them. We will create a single tag whose type is an
423:   // array of two doubles and whose name is "unknowns"
424:   Field default_value = {0.0,0.0};
425:   unknowns_tag_size = sizeof(Field)/sizeof(PetscScalar);
426:   merr = mbint->tag_get_handle("unknowns",unknowns_tag_size,moab::MB_TYPE_DOUBLE,unknowns_tag,
427:                               moab::MB_TAG_DENSE | moab::MB_TAG_CREAT,&default_value);MBERRNM(merr);

429:   // Apply the "unknowns" tag to the vertexes with some initial value...
430:   std::vector<Field> tag_data(my_npts);
431:   for (int i = 0; i < my_npts; i++) {
432:     tag_data[i].u = 1+sin(2*PETSC_PI*vertex_coords[0][i]);
433:     tag_data[i].v = 3;
434:   }

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

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

443:   std::vector<int> global_ids(my_npts);
444:   for (int i = 0; i < my_npts; i++) {
445:     global_ids[i] = i+vstart;
446:   }

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

450:   moab::Range owned_edges;
451:   merr = mbint->get_entities_by_type(0,moab::MBEDGE,owned_edges);MBERRNM(merr);
452:   merr = pcomm->resolve_shared_ents(0,owned_edges,1,0);MBERRNM(merr);

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

458:   // Everything is set up, now just do a tag exchange to update tags
459:   // on all of the shared vertexes:
460:   merr = pcomm->exchange_tags(id_tag,owned_vertexes);MBERRNM(merr);

462:   return(0);
463: }


468: PetscErrorCode create_app_data(_User& user)
469: {
471:   moab::ErrorCode merr;


475:   PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Advection-reaction options",""); {
476:     user.A      = 1;
477:     user.B      = 3;
478:     user.alpha  = 0.02;
479:     user.uleft  = 1;
480:     user.uright = 1;
481:     user.vleft  = 3;
482:     user.vright = 3;
483:     user.npts   = 11;
484:     PetscOptionsReal("-A","Reaction rate","",user.A,&user.A,PETSC_NULL);
485:     PetscOptionsReal("-B","Reaction rate","",user.B,&user.B,PETSC_NULL);
486:     PetscOptionsReal("-alpha","Diffusion coefficient","",user.alpha,&user.alpha,PETSC_NULL);
487:     PetscOptionsReal("-uleft","Dirichlet boundary condition","",user.uleft,&user.uleft,PETSC_NULL);
488:     PetscOptionsReal("-uright","Dirichlet boundary condition","",user.uright,&user.uright,PETSC_NULL);
489:     PetscOptionsReal("-vleft","Dirichlet boundary condition","",user.vleft,&user.vleft,PETSC_NULL);
490:     PetscOptionsReal("-vright","Dirichlet boundary condition","",user.vright,&user.vright,PETSC_NULL);
491:     PetscOptionsInt("-npts","Number of mesh points","",user.npts,&user.npts,PETSC_NULL);
492:   } PetscOptionsEnd();

494:   user.mbint = new moab::Core;
495:   user.pcomm = new moab::ParallelComm(user.mbint,PETSC_COMM_WORLD);

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

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

505:   user.owned_vertexes = new moab::Range;
506:   merr = user.mbint->get_entities_by_type(0,moab::MBVERTEX,*user.owned_vertexes);MBERRNM(merr);
507:   user.all_vertexes = new moab::Range(*user.owned_vertexes);

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

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

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

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

518:   // Do some error checking...make sure that tag_data is in a single
519:   // sequence:
520:   int count;
521:   void *data;
522:   moab::Tag tag;
523:   merr = user.mbint->tag_get_handle("unknowns",tag);MBERRNM(merr);
524:   merr = user.mbint->tag_iterate(tag,user.all_vertexes->begin(),user.all_vertexes->end(),
525:                                  count,data);MBERRNM(merr);
526:   if((unsigned int) count != user.all_vertexes->size()) {
527:     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Tag data not laid out contiguously %i %i",
528:              count,user.all_vertexes->size());
529:   }

531:   return(0);
532: }


537: PetscErrorCode create_matrix(_User &user, DM &dm, Mat *J)
538: {
540:   moab::ErrorCode merr;
541:   moab::Tag ltog_tag;
542:   moab::Range range;

545:   DMMoabGetLocalToGlobalTag(dm,&ltog_tag);
546:   DMMoabGetRange(dm,&range);

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

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

557:   ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, range.size(), gindices,PETSC_COPY_VALUES, &ltog);
558:   MatSetLocalToGlobalMappingBlock(*J,ltog,ltog);

560:   // Clean up:
561:   ISLocalToGlobalMappingDestroy(&ltog);
562:   delete [] gindices;

564:   return(0);
565: }

569: PetscErrorCode destroy_app_data(_User& user)
570: {

573:   delete user.owned_edges;
574:   delete user.owned_vertexes;
575:   delete user.all_vertexes;
576:   delete user.shared_vertexes;
577:   delete user.pcomm;
578:   delete user.mbint;

580:   return(0);
581: }