Actual source code: grid.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: grid.c,v 1.24 2000/07/16 05:40:26 knepley Exp $";
  3: #endif

  5: /*
  6:   This file provides routines for grid vectors (vectors that are associated 
  7:   with grids, possibly with multiple degrees of freedom per node). 

  9:   This component is new; changes in the user interface will be occuring in the
 10:   near future!
 11: */

 13: #include "petscsles.h"            /* For ALE Mesh Support *//*I "petscsles.h" I*/
 14: #include "src/grid/gridimpl.h"    /*I "grid.h" I*//*I "gvec.h" I*/

 16: /* Logging support */
 17: int GRID_COOKIE;
 18: int GRID_Reform, GRID_SetUp;

 20: /*---------------------------------------------- Standard Grid Functions --------------------------------------------*/
 21: #undef  __FUNCT__
 23: /*@
 24:   GridDestroy - Destroys a grid. 

 26:   Collective on Grid

 28:   Input Parameter:
 29: . grid - The grid

 31:   Level: beginner

 33: .keywrods: grid, destroy
 34: .seealso: GridView(), GridCreateGVec(), GridSetUp(), GridRefine()
 35: @*/
 36: int GridDestroy(Grid grid)
 37: {

 42:   if (--grid->refct > 0) return(0);
 43:   (*grid->ops->destroy)(grid);
 44:   PetscLogObjectDestroy(grid);
 45:   PetscHeaderDestroy(grid);
 46:   return(0);
 47: }

 49: #undef  __FUNCT__
 51: /*@
 52:   GridView - Views a grid.

 54:   Collective on Grid

 56:   Input Parameters:
 57: + grid   - The grid
 58: - viewer - The viewer, PETSC_VIEWER_STDOUT_WORLD by default

 60:   Level: beginner

 62: .keywords: grid, view
 63: .seealso: GVecView()
 64: @*/
 65: int GridView(Grid grid, PetscViewer viewer)
 66: {

 71:   /* Do not check return from setup, since we may not have a problem yet */
 72:   if (grid->setupcalled == PETSC_FALSE) GridSetUp(grid);
 73:   if (!viewer) {
 74:     viewer = PETSC_VIEWER_STDOUT_SELF;
 75:   } else {
 77:   }
 78:   (*grid->ops->view)(grid, viewer);
 79:   return(0);
 80: }

 82: EXTERN_C_BEGIN
 83: #undef  __FUNCT__
 85: int GridSerialize_Generic(MPI_Comm comm, Grid *grid, PetscViewer viewer, PetscTruth store) {
 86:   Grid       g;
 87:   Mesh       mesh;
 88:   int        fd, hasParent;
 89:   int        one  = 1;
 90:   int        zero = 0;
 91:   int        field, len;
 92:   char      *type_name = PETSC_NULL;
 93:   PetscTruth match;
 94:   int        ierr;


100:   PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_BINARY, &match);
101:   if (match == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONG, "Must be binary viewer");
102:   PetscViewerBinaryGetDescriptor(viewer, &fd);

104:   if (store) {
105:     g = *grid;
107:     MeshSerialize(comm, &g->mesh, viewer, store);
108:     PetscStrlen(g->type_name, &len);
109:     PetscBinaryWrite(fd, &len,                          1,            PETSC_INT,     0);
110:     PetscBinaryWrite(fd,  g->type_name,                 len,          PETSC_CHAR,    0);

112:     PetscBinaryWrite(fd, &g->numFields,                 1,            PETSC_INT,     0);
113:     for(field = 0; field < g->numFields; field++) {
114:       if (g->fields[field].name != PETSC_NULL) {
115:         PetscStrlen(g->fields[field].name, &len);
116:       } else {
117:         len = 0;
118:       }
119:       PetscBinaryWrite(fd, &len,                        1,            PETSC_INT,     0);
120:       PetscBinaryWrite(fd,  g->fields[field].name,      len,          PETSC_CHAR,    0);
121:       PetscBinaryWrite(fd, &g->fields[field].numComp,   1,            PETSC_INT,     0);
122:       PetscStrlen(g->fields[field].discType, &len);
123:       PetscBinaryWrite(fd, &len,                        1,            PETSC_INT,     0);
124:       PetscBinaryWrite(fd,  g->fields[field].discType,  len,          PETSC_CHAR,    0);
125:       DiscretizationSerialize(g->comm, &g->fields[field].disc, viewer, store);
126:       PetscBinaryWrite(fd, &g->fields[field].isActive,  1,            PETSC_INT,     0);
127:       PetscBinaryWrite(fd, &g->fields[field].isConstrained,      1,   PETSC_INT,     0);
128:       PetscBinaryWrite(fd, &g->fields[field].constraintCompDiff, 1,   PETSC_INT,     0);
129:     }

131:     if (g->gridparent != PETSC_NULL) {
132:       PetscBinaryWrite(fd, &one,                        1,            PETSC_INT,     0);
133:       GridSerialize(comm, &g->gridparent, viewer, store);
134:     } else {
135:       PetscBinaryWrite(fd, &zero,                       1,            PETSC_INT,     0);
136:     }
137:     PetscBinaryWrite(fd, &g->isConstrained,             1,            PETSC_INT,     0);
138: #ifdef NEW_GRID_SERIALIZE
139:     /* This stuff is all made in GridSetUp */
140:     FieldClassSerialize(g->comm, &g->cm, viewer, store);
141:     VarOrderingSerialize(g->cm, &g->order, viewer, store);
142:     FieldClassSerialize(g->comm, &g->constraintCM, viewer, store);
143:     VarOrderingSerialize(g->constraintCM, &g->constraintOrder, viewer, store);
144:     ISSerialize(g->comm, &g->constraintOrdering, viewer, store);
145:     MatSerialize(g->comm, &g->constraintMatrix, viewer, store);
146: #endif

148:     PetscBinaryWrite(fd, &g->reduceAlpha,               1,            PETSC_SCALAR,  0);
149:     PetscBinaryWrite(fd, &g->interpolationType,         1,            PETSC_INT,     0);
150:     PetscBinaryWrite(fd, &g->viewField,                 1,            PETSC_INT,     0);
151:     PetscBinaryWrite(fd, &g->viewComp,                  1,            PETSC_INT,     0);
152:   } else {
153:     MeshSerialize(comm, &mesh, viewer, store);
154:     GridCreate(mesh, &g);
155:     MeshDestroy(mesh);
156:     PetscBinaryRead(fd, &len,                         1,                       PETSC_INT);
157:     if (len) {
158:       PetscMalloc((len+1) * sizeof(char), &type_name);
159:       type_name[len] = 0;
160:     }
161:     PetscBinaryRead(fd,  type_name,                   len,                     PETSC_CHAR);
162:     PetscBinaryRead(fd, &g->numFields,                1,                       PETSC_INT);
163:     g->maxFields = g->numFields;
164:     PetscFree(g->fields);
165:     PetscMalloc(g->numFields * sizeof(Field), &g->fields);
166:     for(field = 0; field < g->numFields; field++) {
167:       PetscBinaryRead(fd, &len,                       1,                       PETSC_INT);
168:       if (len) {
169:         PetscMalloc((len+1) * sizeof(char), &g->fields[field].name);
170:         g->fields[field].name[len] = 0;
171:       }
172:       PetscBinaryRead(fd,  g->fields[field].name,     len,                     PETSC_CHAR);
173:       PetscBinaryRead(fd, &g->fields[field].numComp,  1,                       PETSC_INT);
174:       PetscBinaryRead(fd, &len,                       1,                       PETSC_INT);
175:       PetscMalloc((len+1) * sizeof(DiscretizationType), &g->fields[field].discType);
176:       g->fields[field].discType[len] = 0;
177:       PetscBinaryRead(fd,  g->fields[field].discType, len,                     PETSC_CHAR);
178:       DiscretizationSerialize(g->comm, &g->fields[field].disc, viewer, store);
179:       PetscTypeCompare((PetscObject) g->fields[field].disc, g->fields[field].discType, &match);
180:       if (match == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Incorrect serialization of Discretization");
181:       PetscBinaryRead(fd, &g->fields[field].isActive, 1,                       PETSC_INT);
182:       PetscBinaryRead(fd, &g->fields[field].isConstrained,      1,             PETSC_INT);
183:       PetscBinaryRead(fd, &g->fields[field].constraintCompDiff, 1,             PETSC_INT);
184:     }
185:     GridSetType(g, type_name);
186:     if (len) {
187:       PetscFree(type_name);
188:     }

190:     PetscBinaryRead(fd, &hasParent,                   1,                       PETSC_INT);
191:     if (hasParent) {
192:       GridSerialize(comm, &g->gridparent, viewer, store);
193:     }
194:     PetscBinaryRead(fd, &g->isConstrained,            1,                       PETSC_INT);
195: #ifdef NEW_GRID_SERIALIZE
196:     /* This stuff is all made in GridSetUp */
197:     FieldClassSerialize(g->comm, &g->cm, viewer, store);
198:     VarOrderingSerialize(g->cm, &g->order, viewer, store);
199:     FieldClassSerialize(g->comm, &g->constraintCM, viewer, store);
200:     VarOrderingSerialize(g->constraintCM, &g->constraintOrder, viewer, store);
201:     ISSerialize(g->comm, &g->constraintOrdering, viewer, store);
202:     MatSerialize(g->comm, &g->constraintMatrix, viewer, store);
203: #endif

205:     PetscBinaryRead(fd, &g->reduceAlpha,              1,                       PETSC_SCALAR);
206:     PetscBinaryRead(fd, &g->interpolationType,        1,                       PETSC_INT);
207:     PetscBinaryRead(fd, &g->viewField,                1,                       PETSC_INT);
208:     PetscBinaryRead(fd, &g->viewComp,                 1,                       PETSC_INT);

210:     *grid = g;
211:   }

213:   return(0);
214: }
215: EXTERN_C_END

217: #undef  __FUNCT__
219: /*@
220:   GridSetUp - Set up any required internal data structures for the grid.

222:   Collective on Grid

224:   Input Parameter:
225: . grid - The grid

227:   Level: intermediate

229: .keywords: grid, setup
230: .seealso: GridDestroy(), GridCreateGVec()
231: @*/
232: int GridSetUp(Grid grid)
233: {
234:   LocalVarOrdering reduceLocOrder;
235:   LocalVarOrdering locOrder;
236:   int              sField, tField;
237:   int              bd, bc, op;
238:   int              ierr;

242:   if (grid->setupcalled == PETSC_TRUE) return(0);
243:   grid->setupcalled = PETSC_TRUE;

245:   PetscLogEventBegin(GRID_SetUp, grid, 0, 0, 0);
246:   if (grid->numFields <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "At least one field must be defined on the Grid");
247:   /* Setup boundary sizes */
248:   for(bd = 0; bd < grid->numBd; bd++) {
249:     if (grid->bdSize[bd] != PETSC_NULL) {
250:       PetscFree(grid->bdSize[bd]);
251:     }
252:     PetscMalloc(grid->numFields * sizeof(int), &grid->bdSize[bd]);
253:     PetscLogObjectMemory(grid, grid->numFields * sizeof(int));
254:   }
255:   /* Setup type-specific stuff */
256:   if (grid->ops->gridsetup) {
257:     (*grid->ops->gridsetup)(grid);
258:     if (ierr != 0) PetscFunctionReturn(ierr);
259:   }
260:   /* Setup constraint structures */
261:   (*grid->ops->gridsetupconstraints)(grid, grid->constraintCtx);
262:   /* Setup ghost variable structures */
263:   GridSetupGhostScatter(grid);
264:   /* Setup boundary conditions reductions */
265:   if (grid->reduceSystem == PETSC_TRUE) {
266:     /* Create storage for reduction of Rhs */
267:     if (grid->bdReduceVec != PETSC_NULL) {
268:       GVecDestroy(grid->bdReduceVec);
269:     }
270:     GVecCreateRectangularGhost(grid, grid->reduceOrder, &grid->bdReduceVec);
271:     grid->bdReduceVecCur = grid->bdReduceVec;

273:     /* Hopefully, the user specified the right context by now */
274:     GridCalcBCValues(grid, PETSC_FALSE, grid->reduceContext);

276:     if (grid->reduceElement == PETSC_FALSE) {
277:       /* Create storage for explicit reduction of Rhs */
278:       GVecCreate(grid, &grid->reduceVec);
279:       GMatCreateRectangular(grid, grid->reduceOrder, grid->order, &grid->bdReduceMat);
280:       MatSetOption(grid->bdReduceMat, MAT_NEW_NONZERO_ALLOCATION_ERR);

282:       /* Initialize storage */
283:       MatZeroEntries(grid->bdReduceMat);

285:       /* Evaluate the elimination block */
286:       for(bc = 0; bc < grid->numBC; bc++) {
287:         for(op = 0; op < grid->numMatOps; op++) {
288:           sField = grid->matOps[op].field;
289:           tField = grid->fields[sField].disc->operators[grid->matOps[op].op]->test->field;
290:           if (sField == grid->bc[bc].field) {
291:             LocalVarOrderingCreate(grid, 1, &tField, &locOrder);
292:             LocalVarOrderingCreate(grid, 1, &sField, &reduceLocOrder);
293:             GMatEvaluateOperatorGalerkin(grid->bdReduceMat, PETSC_NULL, 1, &sField, &tField, grid->matOps[op].op,
294:                                                 grid->matOps[op].alpha, MAT_FLUSH_ASSEMBLY, grid->reduceContext);
295: 
296:             LocalVarOrderingDestroy(locOrder);
297:             LocalVarOrderingDestroy(reduceLocOrder);
298:           }
299: #ifdef PETSC_USE_BOPT_g
300:           PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
301: #endif
302:         }
303:       }
304:       for(bc = 0; bc < grid->numPointBC; bc++) {
305:         for(op = 0; op < grid->numMatOps; op++) {
306:           sField = grid->matOps[op].field;
307:           tField = grid->fields[sField].disc->operators[grid->matOps[op].op]->test->field;
308:           if (sField == grid->pointBC[bc].field) {
309:             LocalVarOrderingCreate(grid, 1, &tField, &locOrder);
310:             LocalVarOrderingCreate(grid, 1, &sField, &reduceLocOrder);
311:             GMatEvaluateOperatorGalerkin(grid->bdReduceMat, PETSC_NULL, 1, &sField, &tField, grid->matOps[op].op,
312:                                                 grid->matOps[op].alpha, MAT_FLUSH_ASSEMBLY, grid->reduceContext);
313: 
314:             LocalVarOrderingDestroy(locOrder);
315:             LocalVarOrderingDestroy(reduceLocOrder);
316:           }
317: #ifdef PETSC_USE_BOPT_g
318:           PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
319: #endif
320:         }
321:       }
322:       MatAssemblyBegin(grid->bdReduceMat, MAT_FINAL_ASSEMBLY);
323:       MatAssemblyEnd(grid->bdReduceMat,   MAT_FINAL_ASSEMBLY);
324:     }
325:   }
326:   PetscLogEventEnd(GRID_SetUp, grid, 0, 0, 0);
327:   return(0);
328: }

330: #undef  __FUNCT__
332: /*@
333:   GridSetupGhostScatter - Set up the scatter from a global vector to local values, including ghosts.

335:   Collective on Grid

337:   Input Parameter:
338: . grid - The grid

340:   Level: advanced

342: .keywords: grid, ghost, scatter
343: .seealso: GridDestroy(), GridCreateGVec()
344: @*/
345: int GridSetupGhostScatter(Grid grid)
346: {
347:   PetscTruth isConstrained, explicitConstraints;
348:   int        ierr;

352:   GridSetUp(grid);
353:   /* Destroy old structures */
354:   if (grid->ghostVec     != PETSC_NULL) {
355:     VecDestroy(grid->ghostVec);
356:   }
357:   if (grid->ghostScatter != PETSC_NULL) {
358:     VecScatterDestroy(grid->ghostScatter);
359:   }
360:   /* Setup scatter from global vector to local ghosted vector */
361:   GridIsConstrained(grid, &isConstrained);
362:   GridGetExplicitConstraints(grid, &explicitConstraints);
363:   if ((isConstrained == PETSC_TRUE) && (explicitConstraints == PETSC_TRUE)) {
364:     (*grid->ops->gridsetupghostscatter)(grid, grid->constraintOrder, &grid->ghostVec, &grid->ghostScatter);
365: 
366:   } else {
367:     (*grid->ops->gridsetupghostscatter)(grid, grid->order, &grid->ghostVec, &grid->ghostScatter);
368:   }
369:   return(0);
370: }

372: #undef  __FUNCT__
374: /*@
375:   GridSetupBoundary - Set up any required internal data structures for the boundary.

377:   Collective on Grid

379:   Input Parameter:
380: . grid - The grid

382:   Level: advanced

384: .keywords: grid, setup, boundary
385: .seealso: GridDestroy(), GridCreateGVec()
386: @*/
387: int GridSetupBoundary(Grid grid)
388: {

393:   if (grid->bdSetupCalled == PETSC_TRUE) return(0);
394:   GridSetUp(grid);
395:   if (grid->ops->gridsetupboundary) {
396:     (*grid->ops->gridsetupboundary)(grid);
397:   }
398:   return(0);
399: }

401: #undef  __FUNCT__
403: /*@
404:   GridGetMesh - Returns the context for the underlying mesh.

406:   Not collective

408:   Input Parameter:
409: . grid - The initial grid

411:   Output Parameter:
412: . mesh - The mesh

414:   Level: intermediate

416: .keywords: grid, mesh
417: .seealso: GridDestroy()
418: @*/
419: int GridGetMesh(Grid grid, Mesh *mesh)
420: {
424:   *mesh = grid->mesh;
425:   return(0);
426: }

428: #undef __FUNCT__  
430: /*@
431:   GridGetClassMap - This function returns the default field class map for the grid.

433:   Not collective

435:   Input Parameter:
436: . grid  - The grid

438:   Output Parameter:
439: . order - The default field class map

441:   Level: intermediate

443: .keywords: Grid, get, class map
444: .seealso: GridGetMesh(), GridGetOrder(), FieldClassMapCreate()
445: @*/
446: int GridGetClassMap(Grid grid, FieldClassMap *map)
447: {
451:   if (grid->isConstrained == PETSC_TRUE) {
452:     *map = grid->constraintCM;
453:   } else {
454:     *map = grid->cm;
455:   }
456:   return(0);
457: }

459: #undef __FUNCT__  
461: /*@
462:   GridGetOrder - This function returns the default ordering for the grid.

464:   Not collective

466:   Input Parameter:
467: . grid  - The grid

469:   Output Parameter:
470: . order - The default variable ordering

472:   Level: intermediate

474: .keywords: Grid, get, order
475: .seealso: GridGetLocalOrder(), VarOrderingCreate()
476: @*/
477: int GridGetOrder(Grid grid, VarOrdering *order)
478: {
482:   if (grid->isConstrained == PETSC_TRUE) {
483:     *order = grid->constraintOrder;
484:   } else {
485:     *order = grid->order;
486:   }
487:   return(0);
488: }

490: #undef __FUNCT__  
492: /*@
493:   GridGetLocalOrder - This function returns the local ordering for the grid.

495:   Not collective

497:   Input Parameter:
498: . grid  - The grid

500:   Output Parameter:
501: . order - The local variable ordering

503:   Level: intermediate

505: .keywords: Grid, get, local, order
506: .seealso: GridGetOrder(), LocalVarOrderingCreate()
507: @*/
508: int GridGetLocalOrder(Grid grid, LocalVarOrdering *order)
509: {
513:   *order = grid->locOrder;
514:   return(0);
515: }

517: #undef __FUNCT__  
519: /*@
520:   GridGetDiscretization - This function returns the discretization for the given field.

522:   Not collective

524:   Input Parameters:
525: + grid  - The grid
526: - field - The field

528:   Output Parameter:
529: . disc  - The Discretization

531:   Level: intermediate

533: .keywords: Grid, get, Discretization, field
534: .seealso: GridGetOrder(), DiscretizationCreate()
535: @*/
536: int GridGetDiscretization(Grid grid, int field, Discretization *disc) {
540:   GridValidField(grid, field);
541:   *disc = grid->fields[field].disc;
542:   return(0);
543: }

545: #undef __FUNCT__  
547: /*
548:   GridSetTypeFromOptions - Sets the type of grid from user options.

550:   Collective on Grid

552:   Input Parameter:
553: . grid - The grid

555:   Level: intermediate

557: .keywords: Grid, set, options, database, type
558: .seealso: GridSetFromOptions(), GridSetType()
559: */
560: static int GridSetTypeFromOptions(Grid grid)
561: {
562:   PetscTruth opt;
563:   char      *defaultType;
564:   char       typeName[256];
565:   int        dim = grid->dim;
566:   int        ierr;

569:   if (grid->type_name != PETSC_NULL) {
570:     defaultType = grid->type_name;
571:   } else {
572:     switch(dim) {
573:     case 1:
574:       defaultType = GRID_TRIANGULAR_1D;
575:       break;
576:     case 2:
577:       defaultType = GRID_TRIANGULAR_2D;
578:       break;
579:     default:
580:       SETERRQ1(PETSC_ERR_SUP, "Grid dimension %d is not supported", dim);
581:     }
582:   }

584:   if (GridRegisterAllCalled == PETSC_FALSE) {
585:     GridRegisterAll(PETSC_NULL);
586:   }
587:   PetscOptionsList("-grid_type", "Grid method"," GridSetType", GridList, defaultType, typeName, 256, &opt);
588:   if (opt == PETSC_TRUE) {
589:     GridSetType(grid, typeName);
590:   } else {
591:     GridSetType(grid, defaultType);
592:   }
593:   return(0);
594: }

596: #undef __FUNCT__  
598: /*
599:   GridSetSerializeTypeFromOptions - Sets the type of grid serialization from user options.

601:   Collective on Grid

603:   Input Parameter:
604: . grid - The grid

606:   Level: intermediate

608: .keywords: Grid, set, options, database, type
609: .seealso: GridSetFromOptions(), GridSetSerializeType()
610: */
611: static int GridSetSerializeTypeFromOptions(Grid grid)
612: {
613:   PetscTruth opt;
614:   char      *defaultType;
615:   char       typeName[256];
616:   int        ierr;

619:   if (grid->serialize_name != PETSC_NULL) {
620:     defaultType = grid->serialize_name;
621:   } else {
622:     defaultType = GRID_SER_GENERIC;
623:   }

625:   if (GridSerializeRegisterAllCalled == PETSC_FALSE) {
626:     GridSerializeRegisterAll(PETSC_NULL);
627:   }
628:   PetscOptionsList("-grid_serialize_type", "Grid serialization method"," GridSetSerializeType", GridSerializeList,
629:                           defaultType, typeName, 256, &opt);
630:   if (opt == PETSC_TRUE) {
631:     GridSetSerializeType(grid, typeName);
632:   } else {
633:     GridSetSerializeType(grid, defaultType);
634:   }
635:   return(0);
636: }

638: #undef __FUNCT__  
640: /*@
641:   GridSetFromOptions - Sets various Grid parameters from user options.

643:   Collective on Grid

645:   Input Parameter:
646: . grid - The grid

648:   Notes:  To see all options, run your program with the -help option, or consult the users manual.
649:           Must be called after GridCreate() but before the Grid is used.

651:   Level: intermediate

653: .keywords: Grid, set, options, database
654: .seealso: GridCreate(), GridPrintHelp(), GridSetOptionsPrefix()
655: @*/
656: int GridSetFromOptions(Grid grid)
657: {
658:   Mesh       mesh;
659:   char      *intNames[NUM_INTERPOLATIONS] = {"Local", "Halo", "L2", "L2Local"};
660:   char       intType[256];
661:   PetscTruth opt;
662:   int        ierr;

666:   PetscOptionsBegin(grid->comm, grid->prefix, "Grid options", "Grid");

668:   /* Handle generic grid options */
669:   PetscOptionsHasName(PETSC_NULL, "-help", &opt);
670:   if (opt == PETSC_TRUE) {
671:     GridPrintHelp(grid);
672:   }

674:   /* Constraint options */
675:   PetscOptionsName("-grid_explicit_constraints", "Explicitly form constrained system", "GridGetExplicitConstraints", &grid->explicitConstraints);

677:   /* Interpolation options */
678:   PetscOptionsEList("-grid_int_type", "Interpolation Method", "None", intNames, NUM_INTERPOLATIONS,
679:                            intNames[grid->interpolationType], intType, 256, &opt);
680:   if (opt == PETSC_TRUE) {
681:     PetscTruth islocal, ishalo, isl2, isl2local;

683:     PetscStrcasecmp(intType, intNames[INTERPOLATION_LOCAL],    &islocal);
684:     PetscStrcasecmp(intType, intNames[INTERPOLATION_HALO],     &ishalo);
685:     PetscStrcasecmp(intType, intNames[INTERPOLATION_L2],       &isl2);
686:     PetscStrcasecmp(intType, intNames[INTERPOLATION_L2_LOCAL], &isl2local);
687:     if (islocal == PETSC_TRUE) {
688:       grid->interpolationType = INTERPOLATION_LOCAL;
689:     } else if (ishalo == PETSC_TRUE) {
690:       grid->interpolationType = INTERPOLATION_HALO;
691:     } else if (isl2 == PETSC_TRUE) {
692:       grid->interpolationType = INTERPOLATION_L2;
693:     } else if (isl2local == PETSC_TRUE) {
694:       grid->interpolationType = INTERPOLATION_L2_LOCAL;
695:     }
696:   }

698:   /* Graphics options */
699:   PetscOptionsInt("-grid_view_field", "View Field", "None", grid->viewField, &grid->viewField, &opt);
700:   PetscOptionsInt("-grid_view_comp",  "View Component", "None", grid->viewComp, &grid->viewComp,  &opt);

702:   /* Handle grid type options */
703:   GridSetTypeFromOptions(grid);

705:   /* Handle grid serialization options */
706:   GridSetSerializeTypeFromOptions(grid);

708:   /* Handle specific grid options */
709:   if (grid->ops->setfromoptions != PETSC_NULL) {
710:     (*grid->ops->setfromoptions)(grid);
711:   }

713:   PetscOptionsEnd();

715:   /* Handle subobject options */
716:   GridGetMesh(grid, &mesh);
717:   MeshSetFromOptions(mesh);

719:   GridViewFromOptions(grid, grid->name);
720:   return(0);
721: }

723: #undef  __FUNCT__
725: /*@
726:   GridViewFromOptions - This function visualizes the grid based upon user options.

728:   Collective on Grid

730:   Input Parameter:
731: . grid - The grid

733:   Level: intermediate

735: .keywords: Grid, view, options, database
736: .seealso: GridSetFromOptions(), GridView()
737: @*/
738: int GridViewFromOptions(Grid grid, char *title)
739: {
740:   PetscViewer viewer;
741:   PetscDraw   draw;
742:   PetscTruth  opt;
743:   char       *titleStr;
744:   char        typeName[1024];
745:   char        fileName[PETSC_MAX_PATH_LEN];
746:   int         len;
747:   int         ierr;

750:   PetscOptionsHasName(grid->prefix, "-grid_view", &opt);
751:   if (opt == PETSC_TRUE) {
752:     PetscOptionsGetString(grid->prefix, "-grid_view", typeName, 1024, &opt);
753:     PetscStrlen(typeName, &len);
754:     if (len > 0) {
755:       PetscViewerCreate(grid->comm, &viewer);
756:       PetscViewerSetType(viewer, typeName);
757:       PetscOptionsGetString(grid->prefix, "-grid_view_file", fileName, 1024, &opt);
758:       if (opt == PETSC_TRUE) {
759:         PetscViewerSetFilename(viewer, fileName);
760:       } else {
761:         PetscViewerSetFilename(viewer, grid->name);
762:       }
763:       GridView(grid, viewer);
764:       PetscViewerFlush(viewer);
765:       PetscViewerDestroy(viewer);
766:     } else {
767:       GridView(grid, PETSC_NULL);
768:     }
769:   }
770:   PetscOptionsHasName(grid->prefix, "-grid_view_draw", &opt);
771:   if (opt == PETSC_TRUE) {
772:     PetscViewerDrawOpen(grid->comm, 0, 0, 0, 0, 300, 300, &viewer);
773:     PetscViewerDrawGetDraw(viewer, 0, &draw);
774:     if (title != PETSC_NULL) {
775:       titleStr = title;
776:     } else {
777:       PetscObjectName((PetscObject) grid);                                                         CHKERRQ(ierr) ;
778:       titleStr = grid->name;
779:     }
780:     PetscDrawSetTitle(draw, titleStr);
781:     GridView(grid, viewer);
782:     PetscViewerFlush(viewer);
783:     PetscDrawPause(draw);
784:     PetscViewerDestroy(viewer);
785:   }
786:   return(0);
787: }

789: #undef __FUNCT__  
791: /*@
792:   GridPrintHelp - Prints all options for the grid.

794:   Input Parameter:
795: . grid - The grid

797:   Options Database Keys:
798: $  -help, -h

800:   Level: intermediate

802: .keywords: Grid, help
803: .seealso: GridSetFromOptions()
804: @*/
805: int GridPrintHelp(Grid grid)
806: {
807:   char p[64];
808:   int  ierr;


813:   PetscStrcpy(p, "-");
814:   if (grid->prefix != PETSC_NULL) {
815:     PetscStrcat(p, grid->prefix);
816:   }

818:   (*PetscHelpPrintf)(grid->comm, "Grid options ------------------------------------------------n");
819:   (*PetscHelpPrintf)(grid->comm,"   %sgrid_type <typename> : Sets the grid typen", p);
820:   (*PetscHelpPrintf)(grid->comm,"   %sgrid_view_field <field>: Sets the field to visualizen", p);
821:   (*PetscHelpPrintf)(grid->comm,"   %sgrid_view_comp  <field>: Sets the component to visualizen", p);
822:   return(0);
823: }

825: #undef __FUNCT__  
827: /*@C
828:   GridSetOptionsPrefix - Sets the prefix used for searching for all grid options in the database.

830:   Input Parameters:
831: + grid   - The grid
832: - prefix - The prefix to prepend to all option names

834:   Notes:
835:   A hyphen (-) must NOT be given at the beginning of the prefix name.
836:   The first character of all runtime options is AUTOMATICALLY the hyphen.

838:   Level: intermediate

840: .keywords: grid, set, options, prefix, database
841: .seealso: GridSetFromOptions()
842: @*/
843: int GridSetOptionsPrefix(Grid grid, char *prefix)
844: {

849:   PetscObjectSetOptionsPrefix((PetscObject) grid, prefix);
850:   return(0);
851: }

853: #undef __FUNCT__  
855: /*@C
856:   GridAppendOptionsPrefix - Appends to the prefix used for searching for all grid options in the database.

858:   Collective on Grid

860:   Input Parameters:
861: + grid   - The grid
862: - prefix - The prefix to prepend to all option names

864:   Notes:
865:   A hyphen (-) must NOT be given at the beginning of the prefix name.
866:   The first character of all runtime options is AUTOMATICALLY the hyphen.

868:   Level: intermediate

870: .keywords: grid, append, options, prefix, database
871: .seealso: GridGetOptionsPrefix()
872: @*/
873: int GridAppendOptionsPrefix(Grid grid, char *prefix)
874: {
876: 
879:   PetscObjectAppendOptionsPrefix((PetscObject) grid, prefix);
880:   return(0);
881: }

883: #undef __FUNCT__  
885: /*@
886:   GridGetOptionsPrefix - Sets the prefix used for searching for all grid options in the database.

888:   Input Parameter:
889: . grid   - The grid

891:   Output Parameter:
892: . prefix - A pointer to the prefix string used

894:   Level: intermediate

896: .keywords: grid, get, options, prefix, database
897: .seealso: GridAppendOptionsPrefix()
898: @*/
899: int GridGetOptionsPrefix(Grid grid, char **prefix)
900: {

905:   PetscObjectGetOptionsPrefix((PetscObject) grid, prefix);
906:   return(0);
907: }

909: /*----------------------------------------------- Graphics Functions ------------------------------------------------*/
910: #undef  __FUNCT__
912: /*@C GridSetViewField
913:   This function sets the field component to be used in visualization.

915:   Collective on Grid

917:   Input Parameters:
918: + grid  - The grid
919: . field - The field
920: - comp  - The component

922:   Level: intermediate

924: .keywords grid, view, field
925: .seealso GridGetMeshInfo
926: @*/
927: int GridSetViewField(Grid grid, int field, int comp)
928: {
931:   GridValidField(grid, field);
932:   if ((comp < 0)  || (comp  >= grid->fields[field].numComp)) {
933:     SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid field component %d should be in [0,%d)", comp, grid->fields[field].numComp);
934:   }
935:   grid->viewField = field;
936:   grid->viewComp  = comp;
937:   return(0);
938: }

940: /*------------------------------------------ Interface to Mesh Functions --------------------------------------------*/
941: #undef  __FUNCT__
943: /*@C
944:   GridGetNodeClass - This function returns the class of a given node

946:   Not collective

948:   Input Parameters:
949: + grid   - The grid
950: - node   - The node

952:   Output Parameter:
953: . nclass - The node class

955:   Level: intermediate

957: .keywords grid, mesh, node, class
958: .seealso GridGetMeshInfo
959: @*/
960: int GridGetNodeClass(Grid grid, int node, int *nclass)
961: {
962:   FieldClassMap map;
963:   int           ierr;

968:   GridGetClassMap(grid, &map);
969:   FieldClassMapGetNodeClass(map, node, nclass);
970:   return(0);
971: }

973: #undef  __FUNCT__
975: /*@
976:   GridGetNearestNode - This function returns the node nearest to the given point on which
977:   the field is defined, or -1 if the node is not contained in the local mesh.

979:   Not collective

981:   Input Parameters:
982: + mesh  - The mesh
983: . field - The field
984: - x,y,z - The node coordinates

986:   Output Parameter:
987: . node  - The nearest node

989:   Level: beginner

991: .keywords: grid, node, point location
992: .seealso MeshNearestNode(), MeshLocatePoint()
993: @*/
994: int GridGetNearestNode(Grid grid, int field, double x, double y, double z, int *node)
995: {
996:   FieldClassMap map;
997:   Mesh          mesh;
998:   Partition     p;
999:   PetscTruth    defined;
1000:   double        minDist, dist, nx, ny, nz;
1001:   int           minNode, globalMinNode, allMinNode;
1002:   int           numCorners, startNode, endNode;
1003:   int           elem, corner, nearNode, nclass;
1004:   int           ierr;

1009:   FieldClassMapCreateTriangular2D(grid, 1, &field, &map);
1010:   GridGetMesh(grid, &mesh);
1011:   MeshGetNumCorners(mesh, &numCorners);
1012:   MeshGetPartition(mesh, &p);
1013:   PartitionGetStartNode(p, &startNode);
1014:   PartitionGetEndNode(p, &endNode);
1015:   MeshLocatePoint(mesh, x, y, z, &elem);
1016:   if (elem >= 0) {
1017:     /* Find first node with field defined */
1018:     for(corner = 0, minDist = PETSC_MAX; corner < numCorners; corner++) {
1019:       MeshGetNodeFromElement(mesh, elem, corner, &minNode);
1020:       FieldClassMapGetNodeClass(map, minNode, &nclass);
1021:       FieldClassMapIsDefined(map, field, nclass, &defined);
1022:       if (defined == PETSC_TRUE) {
1023:         MeshGetNodeCoords(mesh, minNode, &nx, &ny, &nz);
1024:         minDist = PetscSqr(MeshPeriodicDiffX(mesh, nx - x)) + PetscSqr(MeshPeriodicDiffY(mesh, ny - y));
1025:         break;
1026:       }
1027:     }
1028:     if (corner == numCorners) SETERRQ2(PETSC_ERR_ARG_WRONG, "Element %d has no node with field %d defined", elem, field);
1029:     /* Find closest node with field defined */
1030:     for(corner = 1; corner < numCorners; corner++) {
1031:       MeshGetNodeFromElement(mesh, elem, corner, &nearNode);
1032:       MeshGetNodeCoords(mesh, nearNode, &nx, &ny, &nz);
1033:       FieldClassMapGetNodeClass(map, nearNode, &nclass);
1034:       FieldClassMapIsDefined(map, field, nclass, &defined);
1035:       dist = PetscSqr(nx - x) + PetscSqr(ny - y);
1036:       if ((defined == PETSC_TRUE) && (dist < minDist)) {
1037:         minDist = dist;
1038:         minNode = nearNode;
1039:       }
1040:     }
1041:     PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode);
1042:   } else {
1043:     minNode       = -1;
1044:     globalMinNode = -1;
1045:   }
1046:   FieldClassMapDestroy(map);

1048:   /* The minimum node might still be a ghost node */
1049:   MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, grid->comm);
1050:   if ((allMinNode >= startNode) && (allMinNode < endNode)) {
1051:     *node = allMinNode - startNode;
1052:   } else {
1053:     *node = -1;
1054:   }
1055:   if (allMinNode < 0) PetscFunctionReturn(1);
1056:   return(0);
1057: }

1059: #undef  __FUNCT__
1061: /*@
1062:   GridGetNearestBdNode - This function returns the boundary node nearest to the given point
1063:   on which the field is defined, or -1 if the node is not contained in the local mesh.

1065:   Not collective

1067:   Input Parameters:
1068: + mesh  - The mesh
1069: . field - The field
1070: - x,y,z - The node coordinates

1072:   Output Parameter:
1073: . node  - The nearest boundary node

1075:   Level: beginner

1077: .keywords: grid, node, point location
1078: .seealso MeshNearestBdNode(), MeshLocatePoint()
1079: @*/
1080: int GridGetNearestBdNode(Grid grid, int field, double x, double y, double z, int *node)
1081: {
1082:   FieldClassMap map;
1083:   Mesh          mesh;
1084:   Partition     p;
1085:   PetscTruth    defined;
1086:   double        minDist, dist, nx, ny, nz;
1087:   int           minNode, globalMinNode, allMinNode;
1088:   int           numCorners, startNode, endNode;
1089:   int           elem, corner, nearNode, nclass, bd;
1090:   int           ierr;

1095:   FieldClassMapCreateTriangular2D(grid, 1, &field, &map);
1096:   GridGetMesh(grid, &mesh);
1097:   MeshGetNumCorners(mesh, &numCorners);
1098:   MeshGetPartition(mesh, &p);
1099:   PartitionGetStartNode(p, &startNode);
1100:   PartitionGetEndNode(p, &endNode);
1101:   MeshLocatePoint(mesh, x, y, z, &elem);
1102:   if (elem >= 0) {
1103:     /* Find first node with field defined */
1104:     for(corner = 0, minDist = PETSC_MAX; corner < numCorners; corner++) {
1105:       MeshGetNodeFromElement(mesh, elem, corner, &minNode);
1106:       MeshGetNodeBoundary(mesh, minNode, &bd);
1107:       FieldClassMapGetNodeClass(map, minNode, &nclass);
1108:       FieldClassMapIsDefined(map, field, nclass, &defined);
1109:       if ((bd != 0) && (defined == PETSC_TRUE)) {
1110:         MeshGetNodeCoords(mesh, minNode, &nx, &ny, &nz);
1111:         minDist = PetscSqr(MeshPeriodicDiffX(mesh, nx - x)) + PetscSqr(MeshPeriodicDiffY(mesh, ny - y));
1112:         break;
1113:       }
1114:     }
1115:     if (corner == numCorners) SETERRQ2(PETSC_ERR_ARG_WRONG, "Element %d has no node with field %d defined", elem, field);
1116:     /* Find closest node with field defined */
1117:     for(corner = 1; corner < numCorners; corner++) {
1118:       MeshGetNodeFromElement(mesh, elem, corner, &nearNode);
1119:       MeshGetNodeCoords(mesh, nearNode, &nx, &ny, &nz);
1120:       MeshGetNodeBoundary(mesh, nearNode, &bd);
1121:       FieldClassMapGetNodeClass(map, nearNode, &nclass);
1122:       FieldClassMapIsDefined(map, field, nclass, &defined);
1123:       dist = PetscSqr(nx - x) + PetscSqr(ny - y);
1124:       if ((bd != 0) && (defined == PETSC_TRUE) && (dist < minDist)) {
1125:         minDist = dist;
1126:         minNode = nearNode;
1127:       }
1128:     }
1129:     PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode);
1130:   } else {
1131:     minNode       = -1;
1132:     globalMinNode = -1;
1133:   }
1134:   FieldClassMapDestroy(map);

1136:   /* The minimum node might still be a ghost node */
1137:   MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, grid->comm);
1138:   if ((allMinNode >= startNode) && (allMinNode < endNode)) {
1139:     *node = allMinNode - startNode;
1140:   } else {
1141:     *node = -1;
1142:   }
1143:   if (allMinNode < 0)
1144:     PetscFunctionReturn(1);
1145:   return(0);
1146: }

1148: #undef  __FUNCT__
1150: /*@
1151:   GridGetBoundarySize - Retrieves the size of the specified boundary for
1152:   a given field.

1154:   Not collective

1156:   Input Parameters:
1157: + grid  - The grid
1158: . bd    - The boundary marker
1159: - field - The field

1161:   Output Parameter:
1162: . size  - The size of the boundary

1164:   Level: intermediate

1166: .keywords: grid, boundary, size
1167: .seealso: GridGetBoundaryNext(), GridGetBoundaryStart()
1168: @*/
1169: int GridGetBoundarySize(Grid grid, int bd, int field, int *size)
1170: {
1171:   int bdIndex;

1177:   GridValidField(grid, field);
1178:   if (grid->setupcalled == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Grid must be setup first")
1179:   ierr  = MeshGetBoundaryIndex(grid->mesh, bd, &bdIndex);
1180:   *size = grid->bdSize[bdIndex][field];
1181:   return(0);
1182: }

1184: #undef  __FUNCT__
1186: /*@
1187:   GridGetBoundaryStart - Retrieves the canonical node number of the first node
1188:   on the specified boundary in a class contained in a given field.

1190:   Not collective

1192:   Input Parameters:
1193: + grid   - The grid
1194: . bd     - The boundary marker
1195: . field  - The field
1196: - ghost  - Flag for including ghost nodes

1198:   Output Parameters:
1199: + node   - The canonical node number
1200: - nclass - The node class

1202:   Level: intermediate

1204: .keywords: grid, boundary, start
1205: .seealso: GridGetBoundaryNext(), MeshGetBoundaryStart
1206: @*/
1207: int GridGetBoundaryStart(Grid grid, int bd, int field, PetscTruth ghost, int *node, int *nclass)
1208: {
1209:   int f;

1216:   if (grid->setupcalled == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Grid must be setup first");
1217:   /* Look for field in default class map */
1218:   for(f = 0; f < grid->cm->numFields; f++) {
1219:     if (grid->cm->fields[f] == field) break;
1220:   }
1221:   if (f == grid->cm->numFields) {
1222:     SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Field %d not present in the default class map", field);
1223:   }
1224:   (*grid->ops->getboundarystart)(grid, bd, f, ghost, grid->cm, node, nclass);
1225:   return(0);
1226: }

1228: #undef  __FUNCT__
1230: /*@
1231:   GridGetBoundaryNext - Retrieves the canonical node number of the next node
1232:   on the specified boundary in a class contained in a given field, with -1
1233:   indicating boundary termination.

1235:   Not collective

1237:   Input Parameters:
1238: + grid   - The grid
1239: . bd     - The boundary marker
1240: . field  - The field
1241: - ghost  - Flag for including ghost nodes

1243:   Output Parameters:
1244: + node   - The canonical node number
1245: - nclass - The node class

1247:   Level: intermediate

1249: .keywords: grid, boundary, next
1250: .seealso: GridGetBoundaryStart(), MeshGetBoundaryNext
1251: @*/
1252: int GridGetBoundaryNext(Grid grid, int bd, int field, PetscTruth ghost, int *node, int *nclass)
1253: {
1254:   int f;

1261:   if (grid->setupcalled == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Grid must be setup first");
1262:   /* Look for field in default class map */
1263:   for(f = 0; f < grid->cm->numFields; f++) {
1264:     if (grid->cm->fields[f] == field) break;
1265:   }
1266:   if (f == grid->cm->numFields) {
1267:     SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Field %d not present in the default class map", field);
1268:   }
1269:   (*grid->ops->getboundarynext)(grid, bd, field, ghost, grid->cm, node, nclass);
1270:   return(0);
1271: }

1273: #undef  __FUNCT__
1275: /*@
1276:   GridRefineMesh - Refine a grid based on area constraints.

1278:   Collective on Grid

1280:   Input Parameters:
1281: + grid    - The initial grid
1282: - area    - A function which gives an area constraint when evaluated inside an element

1284:   Output Parameter:
1285: . newgrid - The refined grid

1287:   Level: advanced

1289: .keywords: grid, refinement
1290: .seealso: MeshRefine(), GridReformMesh()
1291: @*/
1292: int GridRefineMesh(Grid grid, PointFunction area, Grid *newgrid)
1293: {
1294:   Mesh newMesh;
1295:   int  field;
1296:   int  ierr;


1301:   MeshRefine(grid->mesh, area, &newMesh);
1302:   GridCreate(newMesh, newgrid);
1303:   MeshDestroy(newMesh);
1304:   for(field = 0; field < grid->numFields; field++) {
1305:     GridAddField(*newgrid, grid->fields[field].name, grid->fields[field].discType, grid->fields[field].numComp, 0);
1306:   }
1307:   GridSetType(*newgrid, grid->type_name);
1308:   /* Copy registered operators in discretizations */
1309:   /* Set parent: (*newgrid)->gridparent = grid; */

1311:   return(0);
1312: }

1314: /*------------------------------------------ Distributed Data Functions ---------------------------------------------*/
1315: #undef  __FUNCT__
1317: /*@
1318:   GridGlobalToLocalGeneral - Scatters values including ghost variables
1319:   from the given global vector into the given ghost vector.

1321:   Collective on Grid

1323:   Input Parameters:
1324: + grid     - The grid
1325: . vec      - The global vector
1326: . mode     - The insertion mode, either INSERT_VALUES or ADD_VALUES
1327: - scatter  - The scatter

1329:   Output Parameters:
1330: . ghostVec - The local vector

1332:   Level: intermediate

1334: .keywords: grid, scatter, global, local
1335: .seealso: GridGlobalToLocal(), GridLocalToGlobal(), GVecGlobalToLocal()
1336: @*/
1337: int GridGlobalToLocalGeneral(Grid grid, GVec vec, Vec ghostVec, InsertMode mode, VecScatter scatter)
1338: {

1342:   if (mode == ADD_VALUES) SETERRQ(PETSC_ERR_SUP, "No support for adding ghost values");
1343:   VecScatterBegin(vec, ghostVec, mode, SCATTER_FORWARD, scatter);
1344:   VecScatterEnd(vec, ghostVec, mode, SCATTER_FORWARD, scatter);
1345:   return(0);
1346: }

1348: #undef  __FUNCT__
1350: /*@
1351:   GridGlobalToLocal - Scatters values including ghost variables
1352:   from the given global vector into the local grid ghost vector.

1354:   Collective on Grid

1356:   Input Parameters:
1357: + grid - The grid
1358: . mode - The insertion mode, either INSERT_VALUES or ADD_VALUES
1359: - vec  - The grid vector

1361:   Level: intermediate

1363: .keywords: grid, scatter, global, local
1364: .seealso: GridGlobalToLocalGeneral(), GridLocalToGlobal(), GVecGlobalToLocal()
1365: @*/
1366: int GridGlobalToLocal(Grid grid, InsertMode mode, GVec vec)
1367: {
1368:   Grid g;
1369:   int  ierr;


1375:   GVecGetGrid(vec, &g);
1376:   if (grid != g) SETERRQ(PETSC_ERR_ARG_INCOMP, "Vector does not arise from this grid");
1377:   GridGlobalToLocalGeneral(grid, vec, grid->ghostVec, mode, grid->ghostScatter);
1378:   return(0);
1379: }

1381: #undef  __FUNCT__
1383: /*@
1384:   GridLocalToGlobal - Scatters local processor values including ghost variables
1385:   from the local grid ghost vector into the given global vector.

1387:   Collective on Grid

1389:   Input Parameters:
1390: + grid - The grid
1391: . mode - The insertion mode, either SET_VALUES or ADD_VALUES
1392: - vec  - The grid vector

1394:   Level: intermediate

1396: .seealso: GridGlobalToLocal(), GVecLocalToGlobal()
1397: @*/
1398: int GridLocalToGlobal(Grid grid, InsertMode mode, GVec vec)
1399: {
1400:   Grid g;
1401:   int  ierr;


1407:   GVecGetGrid(vec, &g);
1408:   if (grid != g) SETERRQ(PETSC_ERR_ARG_WRONG, "Vector does not arise from this grid");
1409:   VecScatterBegin(grid->ghostVec, vec, mode, SCATTER_REVERSE, grid->ghostScatter);
1410:   VecScatterEnd(grid->ghostVec, vec, mode, SCATTER_REVERSE, grid->ghostScatter);
1411:   return(0);
1412: }

1414: /*---------------------------------------------- InterGrid Functions ------------------------------------------------*/
1415: #undef  __FUNCT__
1417: /*@
1418:   GridInterpolateElementVec - Interpolates a given element vector from one discretization to another.

1420:   Not collective

1422:   Input Parameters:
1423: + grid     - The original grid 
1424: . field    - The original field
1425: . vec      - The original element vector
1426: . newGrid  - The grid defining the new element vector
1427: - newField - The new field

1429:   Output Parameter:
1430: . newVec   - The interpolated element vector

1432:   Level: advanced

1434: .keywords: grid, interpolation
1435: .seealso: ElementVecCreate()
1436: @*/
1437: int GridInterpolateElementVec(Grid grid, int field, ElementVec vec, Grid newGrid, int newField, ElementVec newVec)
1438: {

1446:   GridValidField(grid,    field);
1447:   GridValidField(newGrid, newField);
1448:   if (grid->fields[field].numComp != newGrid->fields[newField].numComp) {
1449:     SETERRQ2(PETSC_ERR_ARG_INCOMP, "Fields have a different number of components %d != %d",
1450:              grid->fields[field].numComp, newGrid->fields[newField].numComp);
1451:   }
1452:   DiscretizationInterpolateElementVec(grid->fields[field].disc, vec, newGrid->fields[newField].disc, newVec);
1453:   return(0);
1454: }

1456: #undef  __FUNCT__
1458: /*@
1459:   GridCreateRestriction - Generates restriction matrix between two grids.

1461:   Input Parameters:
1462: + vf - The fine grid 
1463: - vc - The coarse grid

1465:   Output Parameter:
1466: . gmat - A sparse matrix containing restriction

1468:   Level: advanced

1470: @*/
1471: int GridCreateRestriction(Grid vf, Grid vc, GMat *gmat)
1472: {


1480:   if (vf->setupcalled == PETSC_FALSE) {
1481:     GridSetUp(vf);
1482:   }
1483:   if (vc->setupcalled == PETSC_FALSE) {
1484:     GridSetUp(vc);
1485:   }
1486:   if (!vf->ops->gridcreaterestriction) SETERRQ(PETSC_ERR_SUP, " ");
1487:   (*vf->ops->gridcreaterestriction)(vf, vc, gmat);
1488:   return(0);
1489: }

1491: /*---------------------------------------------- Assembly Functions -------------------------------------------------*/
1492: #undef  __FUNCT__
1494: /*@
1495:   GridLocalToElementGeneral - Scatters values including ghost variables from the given ghost vector
1496:   into the given element vector.

1498:   Not collective

1500:   Input Parameters:
1501: + grid          - The grid
1502: . ghostVec      - The vector of values (including ghsot points)
1503: . reduceVec     - The vector of boundary values
1504: . reduceSystem  - The flag for reducing boundary conditions
1505: . reduceElement - The flag for putting boundary values in the element vector
1506: - vec           - The element vector

1508:   WARNING:
1509:   Make sure that the indices in the element vector are local indices.

1511:   Note:
1512:   If reduceSystem and reduceElement are PTESC_TRUE, then boundary values are
1513:   placed in vec. If reduceElement is PETSC_FALSE, then zero is used instead.

1515:   Level: advanced

1517: .keywords: grid, element, scatter
1518: .seealso: GridLocalToElement(), GridGlobalToLocal()
1519: @*/
1520: int GridLocalToElementGeneral(Grid grid, Vec ghostVec, Vec reduceVec, PetscTruth reduceSystem, PetscTruth reduceElement, ElementVec vec)
1521: {
1522:   PetscScalar *array, *reduceArray;
1523:   PetscScalar  alpha    = -grid->reduceAlpha;
1524:   int          elemSize = vec->reduceSize;
1525:   int          numOverlapVars, numReduceOverlapVars;
1526:   int          var;
1527:   int          ierr;

1530:   VecGetSize(ghostVec, &numOverlapVars);
1531:   VecGetArray(ghostVec, &array);
1532:   if (reduceSystem == PETSC_TRUE) {
1533:     VecGetSize(reduceVec, &numReduceOverlapVars);
1534:     VecGetArray(reduceVec, &reduceArray);
1535:     for(var = 0; var < elemSize; var++) {
1536: #ifdef PETSC_USE_BOPT_g
1537:       if (vec->indices[var] >= numOverlapVars) {
1538:         SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Invalid index %d into ghost vector should be in [0,%d)",
1539:                  vec->indices[var], numOverlapVars);
1540:       }
1541: #endif
1542:       if (vec->indices[var] < 0) {
1543:         if (reduceElement == PETSC_TRUE) {
1544: #ifdef PETSC_USE_BOPT_g
1545:           if (-(vec->indices[var]+1) >= numReduceOverlapVars) {
1546:             SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Invalid index %d into ghost vector should be in [0,%d)",
1547:                      -(vec->indices[var]+1), numReduceOverlapVars);
1548:           }
1549: #endif
1550:           vec->array[var] = alpha*reduceArray[-(vec->indices[var]+1)];
1551:         } else {
1552:           vec->array[var] = 0.0;
1553:         }
1554:       } else {
1555:         vec->array[var] = array[vec->indices[var]];
1556:       }
1557:     }
1558:     VecRestoreArray(reduceVec, &reduceArray);
1559:   } else {
1560:     for(var = 0; var < elemSize; var++) {
1561: #ifdef PETSC_USE_BOPT_g
1562:       if ((vec->indices[var] < 0) ||  (vec->indices[var] >= numOverlapVars)) {
1563:         SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Invalid index %d into ghost vector should be in [0,%d)",
1564:                  vec->indices[var], numOverlapVars);
1565:       }
1566: #endif
1567:       vec->array[var] = array[vec->indices[var]];
1568:     }
1569:   }
1570:   VecRestoreArray(ghostVec, &array);
1571:   return(0);
1572: }

1574: #undef  __FUNCT__
1576: /*@
1577:   GridLocalToElement - Scatters values including ghost variables from the local grid ghost vector
1578:   into the given element vector.

1580:   Not collective

1582:   Input Parameters:
1583: + grid - The grid
1584: - vec  - The element vector

1586:   WARNING:
1587:   Make sure that the indices in the element vector are local indices.

1589:   Level: advanced

1591: .keywords: grid, element, scatter
1592: .seealso: GridLocalToElementGeneral(), GridGlobalToLocal()
1593: @*/
1594: int GridLocalToElement(Grid grid, ElementVec vec)
1595: {

1601:   GridLocalToElementGeneral(grid, grid->ghostVec, grid->bdReduceVecCur, grid->reduceSystem, grid->reduceElement, vec);
1602:   return(0);
1603: }

1605: /*-------------------------------------------- Evaluation Functions -------------------------------------------------*/
1606: #undef  __FUNCT__
1608: /*@C
1609:   GridEvaluateRhs - This function constructs the weak form of the functions and
1610:   operators specified with GridSetRhsFunction(), GridSetRhsOperator(),
1611:   and GridSetRhsNonlinearOperator().

1613:   Collective on Grid

1615:   Input Parameters:
1616: + grid - The grid
1617: . x    - The current iterate
1618: - ctx  - The optional user context

1620:   Output Parameters
1621: . f    - The function value

1623:   Level: beginner

1625: .keywords: grid, rhs
1626: .seealso: GridEvaluateSystemMatrix(), GridSetRhsFunction(), GridSetRhsOperator(), GridSetRhsNonlinearOperator()
1627: @*/
1628: int GridEvaluateRhs(Grid grid, GVec x, GVec f, PetscObject ctx)
1629: {
1630:   Mesh        mesh;
1631:   MeshMover   mover;
1632:   Grid        meshVelGrid;
1633:   PetscScalar alpha = -grid->reduceAlpha;
1634:   int         ierr;

1638:   /* x can be PETSC_NULL */
1640:   /* Setup mesh velocity calculation for ALE variables */
1641:   if (grid->ALEActive == PETSC_TRUE) {
1642:     GridGetMesh(grid, &mesh);
1643:     MeshGetMover(mesh, &mover);
1644:     MeshMoverGetVelocityGrid(mover, &meshVelGrid);
1645:     GridSetUp(meshVelGrid);
1646:   }
1647:   /* Calculate Rhs */
1648:   (*grid->ops->gridevaluaterhs)(grid, x, f, ctx);
1649:   /* Add boundary values */
1650:   if ((grid->reduceSystem == PETSC_TRUE) && (grid->reduceElement == PETSC_FALSE) && (grid->bdReduceMat != PETSC_NULL)) {
1651:     /* Calculate contribution of constrained variables to the Rhs */
1652:     MatMult(grid->bdReduceMat, grid->bdReduceVecCur, grid->reduceVec);
1653:     /* Update Rhs */
1654:     VecAXPY(&alpha, grid->reduceVec, f);
1655:   }
1656:   return(0);
1657: }

1659: #undef  __FUNCT__
1661: /*@C GridEvaluateRhsFunction
1662:   This function constructs the weak form of the functions specified with GridSetRhsFunction().

1664:   Collective on Grid

1666:   Input Parameters:
1667: + grid - The grid
1668: . x    - The current iterate
1669: - ctx  - The optional user context

1671:   Output Parameter:
1672: . f    - The function value

1674:   Level: beginner

1676: .keywords grid, rhs, function
1677: .seealso GridEvaluateRhsOperator(), GridEvaluateRhs(), GridEvaluateSystemMatrix()
1678: @*/
1679: int GridEvaluateRhsFunction(Grid grid, GVec x, GVec f, void *ctx)
1680: {
1681:   PetscTruth oldLinear    = grid->activeOpTypes[1];
1682:   PetscTruth oldNonlinear = grid->activeOpTypes[2];
1683:   int        ierr;

1687:   /* x can be PETSC_NULL */
1689:   /* Turn off computation of operators */
1690:   grid->activeOpTypes[1] = PETSC_FALSE;
1691:   grid->activeOpTypes[2] = PETSC_FALSE;
1692:   GridEvaluateRhs(grid, x, f, (PetscObject) ctx);
1693:   grid->activeOpTypes[1] = oldLinear;
1694:   grid->activeOpTypes[2] = oldNonlinear;
1695:   return(0);
1696: }

1698: #undef  __FUNCT__
1700: /*@C GridEvaluateRhsOperator
1701:   This function constructs the weak form of the operators specified with GridSetRhsOperator().

1703:   Collective on Grid

1705:   Input Parameters:
1706: + grid      - The grid
1707: . x         - The current iterate
1708: . linear    - The flag for including linear operators
1709: . nonlinear - The flag for including nonlinear operators
1710: - ctx       - The optional user context

1712:   Output Parameter:
1713: . f    - The operator value

1715:   Level: beginner

1717: .keywords grid, rhs, operator
1718: .seealso GridEvaluateRhsFunction(), GridEvaluateRhs(), GridEvaluateSystemMatrix()
1719: @*/
1720: int GridEvaluateRhsOperator(Grid grid, GVec x, GVec f, PetscTruth linear, PetscTruth nonlinear, void *ctx)
1721: {
1722:   PetscTruth oldFunction  = grid->activeOpTypes[0];
1723:   PetscTruth oldLinear    = grid->activeOpTypes[1];
1724:   PetscTruth oldNonlinear = grid->activeOpTypes[2];
1725:   int        ierr;

1729:   /* x can be PETSC_NULL */
1731:   /* Turn off computation of operators */
1732:   grid->activeOpTypes[0] = PETSC_FALSE;
1733:   grid->activeOpTypes[1] = linear;
1734:   grid->activeOpTypes[2] = nonlinear;
1735:   GridEvaluateRhs(grid, x, f, (PetscObject) ctx);
1736:   grid->activeOpTypes[0] = oldFunction;
1737:   grid->activeOpTypes[1] = oldLinear;
1738:   grid->activeOpTypes[2] = oldNonlinear;
1739:   return(0);
1740: }

1742: #undef  __FUNCT__
1744: /*@C GridEvaluateSystemMatrix
1745:   This function constructs the weak form of the linear operators
1746:   specified with GridSetMatOperator().

1748:   Collective on Grid

1750:   Input Parameters:
1751: + grid - The grid
1752: . x    - The current iterate
1753: - ctx  - The optional user context

1755:   Output Parameter:
1756: . f    - The function value

1758:   Level: beginner

1760: .keywords: grid, matrix
1761: .seealso: GridEvaluateRhs()
1762: @*/
1763: int GridEvaluateSystemMatrix(Grid grid, GVec x, GMat *J, GMat *M, MatStructure *flag, PetscObject ctx)
1764: {
1765:   Mesh       mesh;
1766:   MeshMover  mover;
1767:   Grid       meshVelGrid;
1768: #ifdef PETSC_HAVE_MATHEMATICA
1769:   PetscTruth opt;
1770: #endif
1771:   int        ierr;

1775:   if (grid->ALEActive == PETSC_TRUE) {
1776:     GridGetMesh(grid, &mesh);
1777:     MeshGetMover(mesh, &mover);
1778:     MeshMoverGetVelocityGrid(mover, &meshVelGrid);
1779:     GridSetUp(meshVelGrid);
1780:   }
1781:   if (grid->isMatrixFree == PETSC_TRUE) {
1782:     /* Should probably set this as the matrix context, but there is no MatShellSetContext() */
1783:     grid->matrixFreeArg     = x;
1784:     grid->matrixFreeContext = ctx;
1785:   } else {
1786:     (*grid->ops->gridevaluatesystemmatrix)(grid, x, J, M, flag, ctx);
1787:   }
1788: #ifdef PETSC_HAVE_MATHEMATICA
1789:   /* Debugging */
1790:   PetscOptionsHasName(PETSC_NULL, "-trace_grid_math", &opt);
1791:   if (opt == PETSC_TRUE) {
1792:     ViewerMathematicaSetName(PETSC_VIEWER_MATHEMATICA_WORLD, "jac");
1793:     MatView(*J, PETSC_VIEWER_MATHEMATICA_WORLD);
1794:     ViewerMathematicaClearName(PETSC_VIEWER_MATHEMATICA_WORLD);
1795:   }
1796: #endif
1797:   return(0);
1798: }

1800: #if 0
1801: /*----------------------------------------- Parallel Communication Functions ----------------------------------------*/
1802: #undef  __FUNCT__
1804: /*@C
1805:   GridGhostExchange - This functions transfers data between local and ghost storage without a predefined mapping.

1807:   Collective on MPI_Comm

1809:   Input Parameters:
1810: + comm         - The communicator
1811: . numGhosts    - The number of ghosts in this domain
1812: . ghostProcs   - The processor from which to obtain each ghost
1813: . ghostIndices - The global index for each ghost
1814: . dataType     - The type of the variables
1815: . firstVar     - The first variable on each processor
1816: . addv         - The insert mode, INSERT_VALUES or ADD_VALUES
1817: - mode         - The direction of the transfer, SCATTER_FORWARD or SCATTER_REVERSE

1819:   Output Paramters:
1820: + locVars      - The local variable array
1821: - ghostVars    - The ghost variables

1823:   Note:
1824:   The data in ghostVars is assumed contiguous and implicitly indexed by the order of
1825:   ghostProcs and ghostIndices. The SCATTER_FORWARD mode will take the requested data
1826:   from locVars and copy it to ghostVars in the order specified by ghostIndices. The
1827:   SCATTER_REVERSE mode will take data from ghostVars and copy it to locVars.

1829:   Level: advanced

1831: .keywords: ghost, exchange, grid
1832: .seealso: GridGlobalToLocal(), GridLocalToGlobal()
1833: @*/
1834: int GridGhostExchange(MPI_Comm comm, int numGhosts, int *ghostProcs, int *ghostIndices, PetscDataType dataType,
1835:                       int *firstVar, InsertMode addv, ScatterMode mode, void *locVars, void *ghostVars)
1836: {
1837:   int         *numSendGhosts; /* The number of ghosts from each domain */
1838:   int         *numRecvGhosts; /* The number of local variables which are ghosts in each domain */
1839:   int         *sumSendGhosts; /* The prefix sums of numSendGhosts */
1840:   int         *sumRecvGhosts; /* The prefix sums of numRecvGhosts */
1841:   int         *offsets;       /* The offset into the send array for each domain */
1842:   int          totSendGhosts; /* The number of ghosts to request variables for */
1843:   int          totRecvGhosts; /* The number of nodes to provide class info about */
1844:   int         *sendIndices;   /* The canonical indices of ghosts in this domain */
1845:   int         *recvIndices;   /* The canonical indices of ghosts to return variables for */
1846:   char        *tempVars;      /* The variables of the requested or submitted ghosts */
1847:   char        *locBytes   = (char *) locVars;
1848:   MPI_Datatype MPIType;
1849:   int          typeSize;
1850: #ifdef PETSC_USE_BOPT_g
1851:   int          numLocVars;
1852: #endif
1853:   int          numProcs, rank;
1854:   int          proc, ghost, locIndex, byte;
1855:   int          ierr;

1858:   /* Initialize communication */
1859:   MPI_Comm_size(comm, &numProcs);
1860:   MPI_Comm_rank(comm, &rank);
1861:   PetscMalloc(numProcs * sizeof(int), &numSendGhosts);
1862:   PetscMalloc(numProcs * sizeof(int), &numRecvGhosts);
1863:   PetscMalloc(numProcs * sizeof(int), &sumSendGhosts);
1864:   PetscMalloc(numProcs * sizeof(int), &sumRecvGhosts);
1865:   PetscMalloc(numProcs * sizeof(int), &offsets);
1866:   PetscMemzero(numSendGhosts,  numProcs * sizeof(int));
1867:   PetscMemzero(numRecvGhosts,  numProcs * sizeof(int));
1868:   PetscMemzero(sumSendGhosts,  numProcs * sizeof(int));
1869:   PetscMemzero(sumRecvGhosts,  numProcs * sizeof(int));
1870:   PetscMemzero(offsets,        numProcs * sizeof(int));
1871: #ifdef PETSC_USE_BOPT_g
1872:   numLocVars = firstVar[rank+1] - firstVar[rank];
1873: #endif

1875:   /* Get number of ghosts needed from each processor */
1876:   for(ghost = 0; ghost < numGhosts; ghost++)
1877:     numSendGhosts[ghostProcs[ghost]]++;

1879:   /* Get number of ghosts to provide variables for */
1880:   MPI_Alltoall(numSendGhosts, 1, MPI_INT, numRecvGhosts, 1, MPI_INT, comm);
1881:   for(proc = 1; proc < numProcs; proc++) {
1882:     sumSendGhosts[proc] = sumSendGhosts[proc-1] + numSendGhosts[proc-1];
1883:     sumRecvGhosts[proc] = sumRecvGhosts[proc-1] + numRecvGhosts[proc-1];
1884:     offsets[proc]       = sumSendGhosts[proc];
1885:   }
1886:   totSendGhosts = sumSendGhosts[numProcs-1] + numSendGhosts[numProcs-1];
1887:   totRecvGhosts = sumRecvGhosts[numProcs-1] + numRecvGhosts[numProcs-1];
1888:   if (numGhosts != totSendGhosts) {
1889:     SETERRQ2(PETSC_ERR_PLIB, "Invalid number of ghosts %d in send, should be %d", totSendGhosts, numGhosts);
1890:   }

1892:   PetscDataTypeGetSize(dataType, &typeSize);
1893:   if (totSendGhosts) {
1894:     sendIndices = (int *)  PetscMalloc(totSendGhosts * sizeof(int)); CHKERRQ(sendIndices);
1895:   }
1896:   if (totRecvGhosts) {
1897:     recvIndices = (int *)  PetscMalloc(totRecvGhosts * sizeof(int)); CHKERRQ(recvIndices);
1898:     tempVars    = (char *) PetscMalloc(totRecvGhosts * typeSize);    CHKERRQ(tempVars);
1899:   }

1901:   /* Must order ghosts by processor */
1902:   for(ghost = 0; ghost < numGhosts; ghost++)
1903:     sendIndices[offsets[ghostProcs[ghost]]++] = ghostIndices[ghost];

1905:   /* Get canonical indices of ghosts to provide variables for */
1906:   MPI_Alltoallv(sendIndices, numSendGhosts, sumSendGhosts, MPI_INT,
1907:                        recvIndices, numRecvGhosts, sumRecvGhosts, MPI_INT, comm);
1908: 

1910:   switch(mode)
1911:   {
1912:   case SCATTER_FORWARD:
1913:     /* Get ghost variables */
1914:     if (addv == INSERT_VALUES) {
1915:       for(ghost = 0; ghost < totRecvGhosts; ghost++)
1916:       {
1917:         locIndex = recvIndices[ghost] - firstVar[rank];
1918: #ifdef PETSC_USE_BOPT_g
1919:         if ((locIndex < 0) || (locIndex >= numLocVars)) {
1920:           SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
1921:         }
1922: #endif
1923:         for(byte = 0; byte < typeSize; byte++)
1924:           tempVars[ghost*typeSize+byte] = locBytes[locIndex*typeSize+byte];
1925:       }
1926:     } else {
1927:       for(ghost = 0; ghost < totRecvGhosts; ghost++)
1928:       {
1929:         locIndex = recvIndices[ghost] - firstVar[rank];
1930: #ifdef PETSC_USE_BOPT_g
1931:         if ((locIndex < 0) || (locIndex >= numLocVars)) {
1932:           SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
1933:         }
1934: #endif
1935:         for(byte = 0; byte < typeSize; byte++)
1936:           tempVars[ghost*typeSize+byte] += locBytes[locIndex*typeSize+byte];
1937:       }
1938:     }

1940:     /* Communicate local variables to ghost storage */
1941:     PetscDataTypeToMPIDataType(dataType, &MPIType);
1942:     MPI_Alltoallv(tempVars,  numRecvGhosts, sumRecvGhosts, MPIType,
1943:                          ghostVars, numSendGhosts, sumSendGhosts, MPIType, comm);
1944: 
1945:     break;
1946:   case SCATTER_REVERSE:
1947:     /* Communicate ghost variables to local storage */
1948:     PetscDataTypeToMPIDataType(dataType, &MPIType);
1949:     MPI_Alltoallv(ghostVars, numSendGhosts, sumSendGhosts, MPIType,
1950:                          tempVars,  numRecvGhosts, sumRecvGhosts, MPIType, comm);
1951: 

1953:     /* Get ghost variables */
1954:     if (addv == INSERT_VALUES) {
1955:       for(ghost = 0; ghost < totRecvGhosts; ghost++)
1956:       {
1957:         locIndex = recvIndices[ghost] - firstVar[rank];
1958: #ifdef PETSC_USE_BOPT_g
1959:         if ((locIndex < 0) || (locIndex >= numLocVars)) {
1960:           SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
1961:         }
1962: #endif
1963:         for(byte = 0; byte < typeSize; byte++)
1964:           locBytes[locIndex*typeSize+byte] = tempVars[ghost*typeSize+byte];
1965:       }
1966:     } else {
1967:       /* There must be a better way to do this -- Ask Bill */
1968:       if (typeSize == sizeof(int)) {
1969:         int *tempInt = (int *) tempVars;
1970:         int *locInt  = (int *) locVars;

1972:         for(ghost = 0; ghost < totRecvGhosts; ghost++) {
1973:           locIndex = recvIndices[ghost] - firstVar[rank];
1974: #ifdef PETSC_USE_BOPT_g
1975:           if ((locIndex < 0) || (locIndex >= numLocVars)) {
1976:             SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
1977:           }
1978: #endif
1979:           locInt[locIndex] += tempInt[ghost];
1980:         }
1981:       } else if (typeSize == sizeof(long int)) {
1982:         long int *tempLongInt = (long int *) tempVars;
1983:         long int *locLongInt  = (long int *) locVars;

1985:         for(ghost = 0; ghost < totRecvGhosts; ghost++) {
1986:           locIndex = recvIndices[ghost] - firstVar[rank];
1987: #ifdef PETSC_USE_BOPT_g
1988:           if ((locIndex < 0) || (locIndex >= numLocVars)) {
1989:             SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
1990:           }
1991: #endif
1992:           locLongInt[locIndex] += tempLongInt[ghost];
1993:         }
1994:       } else {
1995:         for(ghost = 0; ghost < totRecvGhosts; ghost++) {
1996:           locIndex = recvIndices[ghost] - firstVar[rank];
1997: #ifdef PETSC_USE_BOPT_g
1998:           if ((locIndex < 0) || (locIndex >= numLocVars)) {
1999:             SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars);
2000:           }
2001: #endif
2002:           for(byte = 0; byte < typeSize; byte++)
2003:             locBytes[locIndex*typeSize+byte] += tempVars[ghost*typeSize+byte];
2004:         }
2005:       }
2006:     }
2007:     break;
2008:   default:
2009:     SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid scatter mode %d", mode);
2010:   }

2012:   /* Cleanup */
2013:   PetscFree(numSendGhosts);
2014:   PetscFree(numRecvGhosts);
2015:   PetscFree(sumSendGhosts);
2016:   PetscFree(sumRecvGhosts);
2017:   PetscFree(offsets);
2018:   if (totSendGhosts) {
2019:     PetscFree(sendIndices);
2020:   }
2021:   if (totRecvGhosts) {
2022:     PetscFree(recvIndices);
2023:     PetscFree(tempVars);
2024:   }
2025:   return(0);
2026: }
2027: #endif