Actual source code: ex30.cxx
petsc-3.4.5 2014-06-29
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,<og_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, <og);
558: MatSetLocalToGlobalMappingBlock(*J,ltog,ltog);
560: // Clean up:
561: ISLocalToGlobalMappingDestroy(<og);
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: }