#ifdef PETSC_RCS_HEADER static char vcid[] = "$Id: grid.c,v 1.24 2000/07/16 05:40:26 knepley Exp $"; #endif /* This file provides routines for grid vectors (vectors that are associated with grids, possibly with multiple degrees of freedom per node). This component is new; changes in the user interface will be occuring in the near future! */ #include "petscsles.h" /* For ALE Mesh Support *//*I "petscsles.h" I*/ #include "src/grid/gridimpl.h" /*I "grid.h" I*//*I "gvec.h" I*/ /* Logging support */ int GRID_COOKIE; int GRID_Reform, GRID_SetUp; /*---------------------------------------------- Standard Grid Functions --------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridDestroy" /*@ GridDestroy - Destroys a grid. Collective on Grid Input Parameter: . grid - The grid Level: beginner .keywrods: grid, destroy .seealso: GridView(), GridCreateGVec(), GridSetUp(), GridRefine() @*/ int GridDestroy(Grid grid) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); if (--grid->refct > 0) PetscFunctionReturn(0); ierr = (*grid->ops->destroy)(grid); CHKERRQ(ierr); PetscLogObjectDestroy(grid); PetscHeaderDestroy(grid); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridView" /*@ GridView - Views a grid. Collective on Grid Input Parameters: + grid - The grid - viewer - The viewer, PETSC_VIEWER_STDOUT_WORLD by default Level: beginner .keywords: grid, view .seealso: GVecView() @*/ int GridView(Grid grid, PetscViewer viewer) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); /* Do not check return from setup, since we may not have a problem yet */ if (grid->setupcalled == PETSC_FALSE) GridSetUp(grid); if (!viewer) { viewer = PETSC_VIEWER_STDOUT_SELF; } else { PetscValidHeader(viewer); } ierr = (*grid->ops->view)(grid, viewer); CHKERRQ(ierr); PetscFunctionReturn(0); } EXTERN_C_BEGIN #undef __FUNCT__ #define __FUNCT__ "GridSerialize_Generic" int GridSerialize_Generic(MPI_Comm comm, Grid *grid, PetscViewer viewer, PetscTruth store) { Grid g; Mesh mesh; int fd, hasParent; int one = 1; int zero = 0; int field, len; char *type_name = PETSC_NULL; PetscTruth match; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE); PetscValidPointer(grid); ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_BINARY, &match); CHKERRQ(ierr); if (match == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONG, "Must be binary viewer"); ierr = PetscViewerBinaryGetDescriptor(viewer, &fd); CHKERRQ(ierr); if (store) { g = *grid; PetscValidHeaderSpecific(g, GRID_COOKIE); ierr = MeshSerialize(comm, &g->mesh, viewer, store); CHKERRQ(ierr); ierr = PetscStrlen(g->type_name, &len); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &len, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, g->type_name, len, PETSC_CHAR, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &g->numFields, 1, PETSC_INT, 0); CHKERRQ(ierr); for(field = 0; field < g->numFields; field++) { if (g->fields[field].name != PETSC_NULL) { ierr = PetscStrlen(g->fields[field].name, &len); CHKERRQ(ierr); } else { len = 0; } ierr = PetscBinaryWrite(fd, &len, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, g->fields[field].name, len, PETSC_CHAR, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &g->fields[field].numComp, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscStrlen(g->fields[field].discType, &len); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &len, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, g->fields[field].discType, len, PETSC_CHAR, 0); CHKERRQ(ierr); ierr = DiscretizationSerialize(g->comm, &g->fields[field].disc, viewer, store); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &g->fields[field].isActive, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &g->fields[field].isConstrained, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &g->fields[field].constraintCompDiff, 1, PETSC_INT, 0); CHKERRQ(ierr); } if (g->gridparent != PETSC_NULL) { ierr = PetscBinaryWrite(fd, &one, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = GridSerialize(comm, &g->gridparent, viewer, store); CHKERRQ(ierr); } else { ierr = PetscBinaryWrite(fd, &zero, 1, PETSC_INT, 0); CHKERRQ(ierr); } ierr = PetscBinaryWrite(fd, &g->isConstrained, 1, PETSC_INT, 0); CHKERRQ(ierr); #ifdef NEW_GRID_SERIALIZE /* This stuff is all made in GridSetUp */ ierr = FieldClassSerialize(g->comm, &g->cm, viewer, store); CHKERRQ(ierr); ierr = VarOrderingSerialize(g->cm, &g->order, viewer, store); CHKERRQ(ierr); ierr = FieldClassSerialize(g->comm, &g->constraintCM, viewer, store); CHKERRQ(ierr); ierr = VarOrderingSerialize(g->constraintCM, &g->constraintOrder, viewer, store); CHKERRQ(ierr); ierr = ISSerialize(g->comm, &g->constraintOrdering, viewer, store); CHKERRQ(ierr); ierr = MatSerialize(g->comm, &g->constraintMatrix, viewer, store); CHKERRQ(ierr); #endif ierr = PetscBinaryWrite(fd, &g->reduceAlpha, 1, PETSC_SCALAR, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &g->interpolationType, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &g->viewField, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &g->viewComp, 1, PETSC_INT, 0); CHKERRQ(ierr); } else { ierr = MeshSerialize(comm, &mesh, viewer, store); CHKERRQ(ierr); ierr = GridCreate(mesh, &g); ierr = MeshDestroy(mesh); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &len, 1, PETSC_INT); CHKERRQ(ierr); if (len) { ierr = PetscMalloc((len+1) * sizeof(char), &type_name); CHKERRQ(ierr); type_name[len] = 0; } ierr = PetscBinaryRead(fd, type_name, len, PETSC_CHAR); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &g->numFields, 1, PETSC_INT); CHKERRQ(ierr); g->maxFields = g->numFields; ierr = PetscFree(g->fields); CHKERRQ(ierr); ierr = PetscMalloc(g->numFields * sizeof(Field), &g->fields); CHKERRQ(ierr); for(field = 0; field < g->numFields; field++) { ierr = PetscBinaryRead(fd, &len, 1, PETSC_INT); CHKERRQ(ierr); if (len) { ierr = PetscMalloc((len+1) * sizeof(char), &g->fields[field].name); CHKERRQ(ierr); g->fields[field].name[len] = 0; } ierr = PetscBinaryRead(fd, g->fields[field].name, len, PETSC_CHAR); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &g->fields[field].numComp, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &len, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscMalloc((len+1) * sizeof(DiscretizationType), &g->fields[field].discType); CHKERRQ(ierr); g->fields[field].discType[len] = 0; ierr = PetscBinaryRead(fd, g->fields[field].discType, len, PETSC_CHAR); CHKERRQ(ierr); ierr = DiscretizationSerialize(g->comm, &g->fields[field].disc, viewer, store); CHKERRQ(ierr); ierr = PetscTypeCompare((PetscObject) g->fields[field].disc, g->fields[field].discType, &match); CHKERRQ(ierr); if (match == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Incorrect serialization of Discretization"); ierr = PetscBinaryRead(fd, &g->fields[field].isActive, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &g->fields[field].isConstrained, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &g->fields[field].constraintCompDiff, 1, PETSC_INT); CHKERRQ(ierr); } ierr = GridSetType(g, type_name); CHKERRQ(ierr); if (len) { ierr = PetscFree(type_name); CHKERRQ(ierr); } ierr = PetscBinaryRead(fd, &hasParent, 1, PETSC_INT); CHKERRQ(ierr); if (hasParent) { ierr = GridSerialize(comm, &g->gridparent, viewer, store); CHKERRQ(ierr); } ierr = PetscBinaryRead(fd, &g->isConstrained, 1, PETSC_INT); CHKERRQ(ierr); #ifdef NEW_GRID_SERIALIZE /* This stuff is all made in GridSetUp */ ierr = FieldClassSerialize(g->comm, &g->cm, viewer, store); CHKERRQ(ierr); ierr = VarOrderingSerialize(g->cm, &g->order, viewer, store); CHKERRQ(ierr); ierr = FieldClassSerialize(g->comm, &g->constraintCM, viewer, store); CHKERRQ(ierr); ierr = VarOrderingSerialize(g->constraintCM, &g->constraintOrder, viewer, store); CHKERRQ(ierr); ierr = ISSerialize(g->comm, &g->constraintOrdering, viewer, store); CHKERRQ(ierr); ierr = MatSerialize(g->comm, &g->constraintMatrix, viewer, store); CHKERRQ(ierr); #endif ierr = PetscBinaryRead(fd, &g->reduceAlpha, 1, PETSC_SCALAR); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &g->interpolationType, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &g->viewField, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &g->viewComp, 1, PETSC_INT); CHKERRQ(ierr); *grid = g; } PetscFunctionReturn(0); } EXTERN_C_END #undef __FUNCT__ #define __FUNCT__ "GridSetUp" /*@ GridSetUp - Set up any required internal data structures for the grid. Collective on Grid Input Parameter: . grid - The grid Level: intermediate .keywords: grid, setup .seealso: GridDestroy(), GridCreateGVec() @*/ int GridSetUp(Grid grid) { LocalVarOrdering reduceLocOrder; LocalVarOrdering locOrder; int sField, tField; int bd, bc, op; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); if (grid->setupcalled == PETSC_TRUE) PetscFunctionReturn(0); grid->setupcalled = PETSC_TRUE; ierr = PetscLogEventBegin(GRID_SetUp, grid, 0, 0, 0); CHKERRQ(ierr); if (grid->numFields <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "At least one field must be defined on the Grid"); /* Setup boundary sizes */ for(bd = 0; bd < grid->numBd; bd++) { if (grid->bdSize[bd] != PETSC_NULL) { ierr = PetscFree(grid->bdSize[bd]); CHKERRQ(ierr); } ierr = PetscMalloc(grid->numFields * sizeof(int), &grid->bdSize[bd]); CHKERRQ(ierr); PetscLogObjectMemory(grid, grid->numFields * sizeof(int)); } /* Setup type-specific stuff */ if (grid->ops->gridsetup) { ierr = (*grid->ops->gridsetup)(grid); if (ierr != 0) PetscFunctionReturn(ierr); } /* Setup constraint structures */ ierr = (*grid->ops->gridsetupconstraints)(grid, grid->constraintCtx); CHKERRQ(ierr); /* Setup ghost variable structures */ ierr = GridSetupGhostScatter(grid); CHKERRQ(ierr); /* Setup boundary conditions reductions */ if (grid->reduceSystem == PETSC_TRUE) { /* Create storage for reduction of Rhs */ if (grid->bdReduceVec != PETSC_NULL) { ierr = GVecDestroy(grid->bdReduceVec); CHKERRQ(ierr); } ierr = GVecCreateRectangularGhost(grid, grid->reduceOrder, &grid->bdReduceVec); CHKERRQ(ierr); grid->bdReduceVecCur = grid->bdReduceVec; /* Hopefully, the user specified the right context by now */ ierr = GridCalcBCValues(grid, PETSC_FALSE, grid->reduceContext); CHKERRQ(ierr); if (grid->reduceElement == PETSC_FALSE) { /* Create storage for explicit reduction of Rhs */ ierr = GVecCreate(grid, &grid->reduceVec); CHKERRQ(ierr); ierr = GMatCreateRectangular(grid, grid->reduceOrder, grid->order, &grid->bdReduceMat); CHKERRQ(ierr); ierr = MatSetOption(grid->bdReduceMat, MAT_NEW_NONZERO_ALLOCATION_ERR); CHKERRQ(ierr); /* Initialize storage */ ierr = MatZeroEntries(grid->bdReduceMat); CHKERRQ(ierr); /* Evaluate the elimination block */ for(bc = 0; bc < grid->numBC; bc++) { for(op = 0; op < grid->numMatOps; op++) { sField = grid->matOps[op].field; tField = grid->fields[sField].disc->operators[grid->matOps[op].op]->test->field; if (sField == grid->bc[bc].field) { ierr = LocalVarOrderingCreate(grid, 1, &tField, &locOrder); CHKERRQ(ierr); ierr = LocalVarOrderingCreate(grid, 1, &sField, &reduceLocOrder); CHKERRQ(ierr); ierr = GMatEvaluateOperatorGalerkin(grid->bdReduceMat, PETSC_NULL, 1, &sField, &tField, grid->matOps[op].op, grid->matOps[op].alpha, MAT_FLUSH_ASSEMBLY, grid->reduceContext); CHKERRQ(ierr); ierr = LocalVarOrderingDestroy(locOrder); CHKERRQ(ierr); ierr = LocalVarOrderingDestroy(reduceLocOrder); CHKERRQ(ierr); } #ifdef PETSC_USE_BOPT_g ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif } } for(bc = 0; bc < grid->numPointBC; bc++) { for(op = 0; op < grid->numMatOps; op++) { sField = grid->matOps[op].field; tField = grid->fields[sField].disc->operators[grid->matOps[op].op]->test->field; if (sField == grid->pointBC[bc].field) { ierr = LocalVarOrderingCreate(grid, 1, &tField, &locOrder); CHKERRQ(ierr); ierr = LocalVarOrderingCreate(grid, 1, &sField, &reduceLocOrder); CHKERRQ(ierr); ierr = GMatEvaluateOperatorGalerkin(grid->bdReduceMat, PETSC_NULL, 1, &sField, &tField, grid->matOps[op].op, grid->matOps[op].alpha, MAT_FLUSH_ASSEMBLY, grid->reduceContext); CHKERRQ(ierr); ierr = LocalVarOrderingDestroy(locOrder); CHKERRQ(ierr); ierr = LocalVarOrderingDestroy(reduceLocOrder); CHKERRQ(ierr); } #ifdef PETSC_USE_BOPT_g ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif } } ierr = MatAssemblyBegin(grid->bdReduceMat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); ierr = MatAssemblyEnd(grid->bdReduceMat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); } } ierr = PetscLogEventEnd(GRID_SetUp, grid, 0, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetupGhostScatter" /*@ GridSetupGhostScatter - Set up the scatter from a global vector to local values, including ghosts. Collective on Grid Input Parameter: . grid - The grid Level: advanced .keywords: grid, ghost, scatter .seealso: GridDestroy(), GridCreateGVec() @*/ int GridSetupGhostScatter(Grid grid) { PetscTruth isConstrained, explicitConstraints; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); ierr = GridSetUp(grid); CHKERRQ(ierr); /* Destroy old structures */ if (grid->ghostVec != PETSC_NULL) { ierr = VecDestroy(grid->ghostVec); CHKERRQ(ierr); } if (grid->ghostScatter != PETSC_NULL) { ierr = VecScatterDestroy(grid->ghostScatter); CHKERRQ(ierr); } /* Setup scatter from global vector to local ghosted vector */ ierr = GridIsConstrained(grid, &isConstrained); CHKERRQ(ierr); ierr = GridGetExplicitConstraints(grid, &explicitConstraints); CHKERRQ(ierr); if ((isConstrained == PETSC_TRUE) && (explicitConstraints == PETSC_TRUE)) { ierr = (*grid->ops->gridsetupghostscatter)(grid, grid->constraintOrder, &grid->ghostVec, &grid->ghostScatter); CHKERRQ(ierr); } else { ierr = (*grid->ops->gridsetupghostscatter)(grid, grid->order, &grid->ghostVec, &grid->ghostScatter); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetupBoundary" /*@ GridSetupBoundary - Set up any required internal data structures for the boundary. Collective on Grid Input Parameter: . grid - The grid Level: advanced .keywords: grid, setup, boundary .seealso: GridDestroy(), GridCreateGVec() @*/ int GridSetupBoundary(Grid grid) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); if (grid->bdSetupCalled == PETSC_TRUE) PetscFunctionReturn(0); ierr = GridSetUp(grid); CHKERRQ(ierr); if (grid->ops->gridsetupboundary) { ierr = (*grid->ops->gridsetupboundary)(grid); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetMesh" /*@ GridGetMesh - Returns the context for the underlying mesh. Not collective Input Parameter: . grid - The initial grid Output Parameter: . mesh - The mesh Level: intermediate .keywords: grid, mesh .seealso: GridDestroy() @*/ int GridGetMesh(Grid grid, Mesh *mesh) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(mesh); *mesh = grid->mesh; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetClassMap" /*@ GridGetClassMap - This function returns the default field class map for the grid. Not collective Input Parameter: . grid - The grid Output Parameter: . order - The default field class map Level: intermediate .keywords: Grid, get, class map .seealso: GridGetMesh(), GridGetOrder(), FieldClassMapCreate() @*/ int GridGetClassMap(Grid grid, FieldClassMap *map) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(map); if (grid->isConstrained == PETSC_TRUE) { *map = grid->constraintCM; } else { *map = grid->cm; } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetOrder" /*@ GridGetOrder - This function returns the default ordering for the grid. Not collective Input Parameter: . grid - The grid Output Parameter: . order - The default variable ordering Level: intermediate .keywords: Grid, get, order .seealso: GridGetLocalOrder(), VarOrderingCreate() @*/ int GridGetOrder(Grid grid, VarOrdering *order) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(order); if (grid->isConstrained == PETSC_TRUE) { *order = grid->constraintOrder; } else { *order = grid->order; } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetLocalOrder" /*@ GridGetLocalOrder - This function returns the local ordering for the grid. Not collective Input Parameter: . grid - The grid Output Parameter: . order - The local variable ordering Level: intermediate .keywords: Grid, get, local, order .seealso: GridGetOrder(), LocalVarOrderingCreate() @*/ int GridGetLocalOrder(Grid grid, LocalVarOrdering *order) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(order); *order = grid->locOrder; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetDiscretization" /*@ GridGetDiscretization - This function returns the discretization for the given field. Not collective Input Parameters: + grid - The grid - field - The field Output Parameter: . disc - The Discretization Level: intermediate .keywords: Grid, get, Discretization, field .seealso: GridGetOrder(), DiscretizationCreate() @*/ int GridGetDiscretization(Grid grid, int field, Discretization *disc) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(disc); GridValidField(grid, field); *disc = grid->fields[field].disc; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetTypeFromOptions" /* GridSetTypeFromOptions - Sets the type of grid from user options. Collective on Grid Input Parameter: . grid - The grid Level: intermediate .keywords: Grid, set, options, database, type .seealso: GridSetFromOptions(), GridSetType() */ static int GridSetTypeFromOptions(Grid grid) { PetscTruth opt; char *defaultType; char typeName[256]; int dim = grid->dim; int ierr; PetscFunctionBegin; if (grid->type_name != PETSC_NULL) { defaultType = grid->type_name; } else { switch(dim) { case 1: defaultType = GRID_TRIANGULAR_1D; break; case 2: defaultType = GRID_TRIANGULAR_2D; break; default: SETERRQ1(PETSC_ERR_SUP, "Grid dimension %d is not supported", dim); } } if (GridRegisterAllCalled == PETSC_FALSE) { ierr = GridRegisterAll(PETSC_NULL); CHKERRQ(ierr); } ierr = PetscOptionsList("-grid_type", "Grid method"," GridSetType", GridList, defaultType, typeName, 256, &opt);CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = GridSetType(grid, typeName); CHKERRQ(ierr); } else { ierr = GridSetType(grid, defaultType); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetSerializeTypeFromOptions" /* GridSetSerializeTypeFromOptions - Sets the type of grid serialization from user options. Collective on Grid Input Parameter: . grid - The grid Level: intermediate .keywords: Grid, set, options, database, type .seealso: GridSetFromOptions(), GridSetSerializeType() */ static int GridSetSerializeTypeFromOptions(Grid grid) { PetscTruth opt; char *defaultType; char typeName[256]; int ierr; PetscFunctionBegin; if (grid->serialize_name != PETSC_NULL) { defaultType = grid->serialize_name; } else { defaultType = GRID_SER_GENERIC; } if (GridSerializeRegisterAllCalled == PETSC_FALSE) { ierr = GridSerializeRegisterAll(PETSC_NULL); CHKERRQ(ierr); } ierr = PetscOptionsList("-grid_serialize_type", "Grid serialization method"," GridSetSerializeType", GridSerializeList, defaultType, typeName, 256, &opt);CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = GridSetSerializeType(grid, typeName); CHKERRQ(ierr); } else { ierr = GridSetSerializeType(grid, defaultType); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetFromOptions" /*@ GridSetFromOptions - Sets various Grid parameters from user options. Collective on Grid Input Parameter: . grid - The grid Notes: To see all options, run your program with the -help option, or consult the users manual. Must be called after GridCreate() but before the Grid is used. Level: intermediate .keywords: Grid, set, options, database .seealso: GridCreate(), GridPrintHelp(), GridSetOptionsPrefix() @*/ int GridSetFromOptions(Grid grid) { Mesh mesh; char *intNames[NUM_INTERPOLATIONS] = {"Local", "Halo", "L2", "L2Local"}; char intType[256]; PetscTruth opt; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); ierr = PetscOptionsBegin(grid->comm, grid->prefix, "Grid options", "Grid"); CHKERRQ(ierr); /* Handle generic grid options */ ierr = PetscOptionsHasName(PETSC_NULL, "-help", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = GridPrintHelp(grid); CHKERRQ(ierr); } /* Constraint options */ ierr = PetscOptionsName("-grid_explicit_constraints", "Explicitly form constrained system", "GridGetExplicitConstraints", &grid->explicitConstraints);CHKERRQ(ierr); /* Interpolation options */ ierr = PetscOptionsEList("-grid_int_type", "Interpolation Method", "None", intNames, NUM_INTERPOLATIONS, intNames[grid->interpolationType], intType, 256, &opt);CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscTruth islocal, ishalo, isl2, isl2local; ierr = PetscStrcasecmp(intType, intNames[INTERPOLATION_LOCAL], &islocal); CHKERRQ(ierr); ierr = PetscStrcasecmp(intType, intNames[INTERPOLATION_HALO], &ishalo); CHKERRQ(ierr); ierr = PetscStrcasecmp(intType, intNames[INTERPOLATION_L2], &isl2); CHKERRQ(ierr); ierr = PetscStrcasecmp(intType, intNames[INTERPOLATION_L2_LOCAL], &isl2local); CHKERRQ(ierr); if (islocal == PETSC_TRUE) { grid->interpolationType = INTERPOLATION_LOCAL; } else if (ishalo == PETSC_TRUE) { grid->interpolationType = INTERPOLATION_HALO; } else if (isl2 == PETSC_TRUE) { grid->interpolationType = INTERPOLATION_L2; } else if (isl2local == PETSC_TRUE) { grid->interpolationType = INTERPOLATION_L2_LOCAL; } } /* Graphics options */ ierr = PetscOptionsInt("-grid_view_field", "View Field", "None", grid->viewField, &grid->viewField, &opt);CHKERRQ(ierr); ierr = PetscOptionsInt("-grid_view_comp", "View Component", "None", grid->viewComp, &grid->viewComp, &opt);CHKERRQ(ierr); /* Handle grid type options */ ierr = GridSetTypeFromOptions(grid); CHKERRQ(ierr); /* Handle grid serialization options */ ierr = GridSetSerializeTypeFromOptions(grid); CHKERRQ(ierr); /* Handle specific grid options */ if (grid->ops->setfromoptions != PETSC_NULL) { ierr = (*grid->ops->setfromoptions)(grid); CHKERRQ(ierr); } ierr = PetscOptionsEnd(); CHKERRQ(ierr); /* Handle subobject options */ ierr = GridGetMesh(grid, &mesh); CHKERRQ(ierr); ierr = MeshSetFromOptions(mesh); CHKERRQ(ierr); ierr = GridViewFromOptions(grid, grid->name); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridViewFromOptions" /*@ GridViewFromOptions - This function visualizes the grid based upon user options. Collective on Grid Input Parameter: . grid - The grid Level: intermediate .keywords: Grid, view, options, database .seealso: GridSetFromOptions(), GridView() @*/ int GridViewFromOptions(Grid grid, char *title) { PetscViewer viewer; PetscDraw draw; PetscTruth opt; char *titleStr; char typeName[1024]; char fileName[PETSC_MAX_PATH_LEN]; int len; int ierr; PetscFunctionBegin; ierr = PetscOptionsHasName(grid->prefix, "-grid_view", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = PetscOptionsGetString(grid->prefix, "-grid_view", typeName, 1024, &opt); CHKERRQ(ierr); ierr = PetscStrlen(typeName, &len); CHKERRQ(ierr); if (len > 0) { ierr = PetscViewerCreate(grid->comm, &viewer); CHKERRQ(ierr); ierr = PetscViewerSetType(viewer, typeName); CHKERRQ(ierr); ierr = PetscOptionsGetString(grid->prefix, "-grid_view_file", fileName, 1024, &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = PetscViewerSetFilename(viewer, fileName); CHKERRQ(ierr); } else { ierr = PetscViewerSetFilename(viewer, grid->name); CHKERRQ(ierr); } ierr = GridView(grid, viewer); CHKERRQ(ierr); ierr = PetscViewerFlush(viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(viewer); CHKERRQ(ierr); } else { ierr = GridView(grid, PETSC_NULL); CHKERRQ(ierr); } } ierr = PetscOptionsHasName(grid->prefix, "-grid_view_draw", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = PetscViewerDrawOpen(grid->comm, 0, 0, 0, 0, 300, 300, &viewer); CHKERRQ(ierr); ierr = PetscViewerDrawGetDraw(viewer, 0, &draw); CHKERRQ(ierr); if (title != PETSC_NULL) { titleStr = title; } else { ierr = PetscObjectName((PetscObject) grid); CHKERRQ(ierr) ; titleStr = grid->name; } ierr = PetscDrawSetTitle(draw, titleStr); CHKERRQ(ierr); ierr = GridView(grid, viewer); CHKERRQ(ierr); ierr = PetscViewerFlush(viewer); CHKERRQ(ierr); ierr = PetscDrawPause(draw); CHKERRQ(ierr); ierr = PetscViewerDestroy(viewer); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridPrintHelp" /*@ GridPrintHelp - Prints all options for the grid. Input Parameter: . grid - The grid Options Database Keys: $ -help, -h Level: intermediate .keywords: Grid, help .seealso: GridSetFromOptions() @*/ int GridPrintHelp(Grid grid) { char p[64]; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); ierr = PetscStrcpy(p, "-"); CHKERRQ(ierr); if (grid->prefix != PETSC_NULL) { ierr = PetscStrcat(p, grid->prefix); CHKERRQ(ierr); } (*PetscHelpPrintf)(grid->comm, "Grid options ------------------------------------------------\n"); (*PetscHelpPrintf)(grid->comm," %sgrid_type : Sets the grid type\n", p); (*PetscHelpPrintf)(grid->comm," %sgrid_view_field : Sets the field to visualize\n", p); (*PetscHelpPrintf)(grid->comm," %sgrid_view_comp : Sets the component to visualize\n", p); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetOptionsPrefix" /*@C GridSetOptionsPrefix - Sets the prefix used for searching for all grid options in the database. Input Parameters: + grid - The grid - prefix - The prefix to prepend to all option names Notes: A hyphen (-) must NOT be given at the beginning of the prefix name. The first character of all runtime options is AUTOMATICALLY the hyphen. Level: intermediate .keywords: grid, set, options, prefix, database .seealso: GridSetFromOptions() @*/ int GridSetOptionsPrefix(Grid grid, char *prefix) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); ierr = PetscObjectSetOptionsPrefix((PetscObject) grid, prefix); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridAppendOptionsPrefix" /*@C GridAppendOptionsPrefix - Appends to the prefix used for searching for all grid options in the database. Collective on Grid Input Parameters: + grid - The grid - prefix - The prefix to prepend to all option names Notes: A hyphen (-) must NOT be given at the beginning of the prefix name. The first character of all runtime options is AUTOMATICALLY the hyphen. Level: intermediate .keywords: grid, append, options, prefix, database .seealso: GridGetOptionsPrefix() @*/ int GridAppendOptionsPrefix(Grid grid, char *prefix) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); ierr = PetscObjectAppendOptionsPrefix((PetscObject) grid, prefix); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetOptionsPrefix" /*@ GridGetOptionsPrefix - Sets the prefix used for searching for all grid options in the database. Input Parameter: . grid - The grid Output Parameter: . prefix - A pointer to the prefix string used Level: intermediate .keywords: grid, get, options, prefix, database .seealso: GridAppendOptionsPrefix() @*/ int GridGetOptionsPrefix(Grid grid, char **prefix) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); ierr = PetscObjectGetOptionsPrefix((PetscObject) grid, prefix); CHKERRQ(ierr); PetscFunctionReturn(0); } /*----------------------------------------------- Graphics Functions ------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridSetViewField" /*@C GridSetViewField This function sets the field component to be used in visualization. Collective on Grid Input Parameters: + grid - The grid . field - The field - comp - The component Level: intermediate .keywords grid, view, field .seealso GridGetMeshInfo @*/ int GridSetViewField(Grid grid, int field, int comp) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); GridValidField(grid, field); if ((comp < 0) || (comp >= grid->fields[field].numComp)) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid field component %d should be in [0,%d)", comp, grid->fields[field].numComp); } grid->viewField = field; grid->viewComp = comp; PetscFunctionReturn(0); } /*------------------------------------------ Interface to Mesh Functions --------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridGetNodeClass" /*@C GridGetNodeClass - This function returns the class of a given node Not collective Input Parameters: + grid - The grid - node - The node Output Parameter: . nclass - The node class Level: intermediate .keywords grid, mesh, node, class .seealso GridGetMeshInfo @*/ int GridGetNodeClass(Grid grid, int node, int *nclass) { FieldClassMap map; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidIntPointer(nclass); ierr = GridGetClassMap(grid, &map); CHKERRQ(ierr); ierr = FieldClassMapGetNodeClass(map, node, nclass); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetNearestNode" /*@ GridGetNearestNode - This function returns the node nearest to the given point on which the field is defined, or -1 if the node is not contained in the local mesh. Not collective Input Parameters: + mesh - The mesh . field - The field - x,y,z - The node coordinates Output Parameter: . node - The nearest node Level: beginner .keywords: grid, node, point location .seealso MeshNearestNode(), MeshLocatePoint() @*/ int GridGetNearestNode(Grid grid, int field, double x, double y, double z, int *node) { FieldClassMap map; Mesh mesh; Partition p; PetscTruth defined; double minDist, dist, nx, ny, nz; int minNode, globalMinNode, allMinNode; int numCorners, startNode, endNode; int elem, corner, nearNode, nclass; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidIntPointer(node); ierr = FieldClassMapCreateTriangular2D(grid, 1, &field, &map); CHKERRQ(ierr); ierr = GridGetMesh(grid, &mesh); CHKERRQ(ierr); ierr = MeshGetNumCorners(mesh, &numCorners); CHKERRQ(ierr); ierr = MeshGetPartition(mesh, &p); CHKERRQ(ierr); ierr = PartitionGetStartNode(p, &startNode); CHKERRQ(ierr); ierr = PartitionGetEndNode(p, &endNode); CHKERRQ(ierr); ierr = MeshLocatePoint(mesh, x, y, z, &elem); CHKERRQ(ierr); if (elem >= 0) { /* Find first node with field defined */ for(corner = 0, minDist = PETSC_MAX; corner < numCorners; corner++) { ierr = MeshGetNodeFromElement(mesh, elem, corner, &minNode); CHKERRQ(ierr); ierr = FieldClassMapGetNodeClass(map, minNode, &nclass); CHKERRQ(ierr); ierr = FieldClassMapIsDefined(map, field, nclass, &defined); CHKERRQ(ierr); if (defined == PETSC_TRUE) { ierr = MeshGetNodeCoords(mesh, minNode, &nx, &ny, &nz); CHKERRQ(ierr); minDist = PetscSqr(MeshPeriodicDiffX(mesh, nx - x)) + PetscSqr(MeshPeriodicDiffY(mesh, ny - y)); break; } } if (corner == numCorners) SETERRQ2(PETSC_ERR_ARG_WRONG, "Element %d has no node with field %d defined", elem, field); /* Find closest node with field defined */ for(corner = 1; corner < numCorners; corner++) { ierr = MeshGetNodeFromElement(mesh, elem, corner, &nearNode); CHKERRQ(ierr); ierr = MeshGetNodeCoords(mesh, nearNode, &nx, &ny, &nz); CHKERRQ(ierr); ierr = FieldClassMapGetNodeClass(map, nearNode, &nclass); CHKERRQ(ierr); ierr = FieldClassMapIsDefined(map, field, nclass, &defined); CHKERRQ(ierr); dist = PetscSqr(nx - x) + PetscSqr(ny - y); if ((defined == PETSC_TRUE) && (dist < minDist)) { minDist = dist; minNode = nearNode; } } ierr = PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode); CHKERRQ(ierr); } else { minNode = -1; globalMinNode = -1; } ierr = FieldClassMapDestroy(map); CHKERRQ(ierr); /* The minimum node might still be a ghost node */ ierr = MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, grid->comm); CHKERRQ(ierr); if ((allMinNode >= startNode) && (allMinNode < endNode)) { *node = allMinNode - startNode; } else { *node = -1; } if (allMinNode < 0) PetscFunctionReturn(1); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetNearestBdNode" /*@ GridGetNearestBdNode - This function returns the boundary node nearest to the given point on which the field is defined, or -1 if the node is not contained in the local mesh. Not collective Input Parameters: + mesh - The mesh . field - The field - x,y,z - The node coordinates Output Parameter: . node - The nearest boundary node Level: beginner .keywords: grid, node, point location .seealso MeshNearestBdNode(), MeshLocatePoint() @*/ int GridGetNearestBdNode(Grid grid, int field, double x, double y, double z, int *node) { FieldClassMap map; Mesh mesh; Partition p; PetscTruth defined; double minDist, dist, nx, ny, nz; int minNode, globalMinNode, allMinNode; int numCorners, startNode, endNode; int elem, corner, nearNode, nclass, bd; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidIntPointer(node); ierr = FieldClassMapCreateTriangular2D(grid, 1, &field, &map); CHKERRQ(ierr); ierr = GridGetMesh(grid, &mesh); CHKERRQ(ierr); ierr = MeshGetNumCorners(mesh, &numCorners); CHKERRQ(ierr); ierr = MeshGetPartition(mesh, &p); CHKERRQ(ierr); ierr = PartitionGetStartNode(p, &startNode); CHKERRQ(ierr); ierr = PartitionGetEndNode(p, &endNode); CHKERRQ(ierr); ierr = MeshLocatePoint(mesh, x, y, z, &elem); CHKERRQ(ierr); if (elem >= 0) { /* Find first node with field defined */ for(corner = 0, minDist = PETSC_MAX; corner < numCorners; corner++) { ierr = MeshGetNodeFromElement(mesh, elem, corner, &minNode); CHKERRQ(ierr); ierr = MeshGetNodeBoundary(mesh, minNode, &bd); CHKERRQ(ierr); ierr = FieldClassMapGetNodeClass(map, minNode, &nclass); CHKERRQ(ierr); ierr = FieldClassMapIsDefined(map, field, nclass, &defined); CHKERRQ(ierr); if ((bd != 0) && (defined == PETSC_TRUE)) { ierr = MeshGetNodeCoords(mesh, minNode, &nx, &ny, &nz); CHKERRQ(ierr); minDist = PetscSqr(MeshPeriodicDiffX(mesh, nx - x)) + PetscSqr(MeshPeriodicDiffY(mesh, ny - y)); break; } } if (corner == numCorners) SETERRQ2(PETSC_ERR_ARG_WRONG, "Element %d has no node with field %d defined", elem, field); /* Find closest node with field defined */ for(corner = 1; corner < numCorners; corner++) { ierr = MeshGetNodeFromElement(mesh, elem, corner, &nearNode); CHKERRQ(ierr); ierr = MeshGetNodeCoords(mesh, nearNode, &nx, &ny, &nz); CHKERRQ(ierr); ierr = MeshGetNodeBoundary(mesh, nearNode, &bd); CHKERRQ(ierr); ierr = FieldClassMapGetNodeClass(map, nearNode, &nclass); CHKERRQ(ierr); ierr = FieldClassMapIsDefined(map, field, nclass, &defined); CHKERRQ(ierr); dist = PetscSqr(nx - x) + PetscSqr(ny - y); if ((bd != 0) && (defined == PETSC_TRUE) && (dist < minDist)) { minDist = dist; minNode = nearNode; } } ierr = PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode); CHKERRQ(ierr); } else { minNode = -1; globalMinNode = -1; } ierr = FieldClassMapDestroy(map); CHKERRQ(ierr); /* The minimum node might still be a ghost node */ ierr = MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, grid->comm); CHKERRQ(ierr); if ((allMinNode >= startNode) && (allMinNode < endNode)) { *node = allMinNode - startNode; } else { *node = -1; } if (allMinNode < 0) PetscFunctionReturn(1); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetBoundarySize" /*@ GridGetBoundarySize - Retrieves the size of the specified boundary for a given field. Not collective Input Parameters: + grid - The grid . bd - The boundary marker - field - The field Output Parameter: . size - The size of the boundary Level: intermediate .keywords: grid, boundary, size .seealso: GridGetBoundaryNext(), GridGetBoundaryStart() @*/ int GridGetBoundarySize(Grid grid, int bd, int field, int *size) { int bdIndex; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidIntPointer(size); GridValidField(grid, field); if (grid->setupcalled == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Grid must be setup first") ierr = MeshGetBoundaryIndex(grid->mesh, bd, &bdIndex); CHKERRQ(ierr); *size = grid->bdSize[bdIndex][field]; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetBoundaryStart" /*@ GridGetBoundaryStart - Retrieves the canonical node number of the first node on the specified boundary in a class contained in a given field. Not collective Input Parameters: + grid - The grid . bd - The boundary marker . field - The field - ghost - Flag for including ghost nodes Output Parameters: + node - The canonical node number - nclass - The node class Level: intermediate .keywords: grid, boundary, start .seealso: GridGetBoundaryNext(), MeshGetBoundaryStart @*/ int GridGetBoundaryStart(Grid grid, int bd, int field, PetscTruth ghost, int *node, int *nclass) { int f; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidIntPointer(node); PetscValidIntPointer(nclass); if (grid->setupcalled == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Grid must be setup first"); /* Look for field in default class map */ for(f = 0; f < grid->cm->numFields; f++) { if (grid->cm->fields[f] == field) break; } if (f == grid->cm->numFields) { SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Field %d not present in the default class map", field); } ierr = (*grid->ops->getboundarystart)(grid, bd, f, ghost, grid->cm, node, nclass); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetBoundaryNext" /*@ GridGetBoundaryNext - Retrieves the canonical node number of the next node on the specified boundary in a class contained in a given field, with -1 indicating boundary termination. Not collective Input Parameters: + grid - The grid . bd - The boundary marker . field - The field - ghost - Flag for including ghost nodes Output Parameters: + node - The canonical node number - nclass - The node class Level: intermediate .keywords: grid, boundary, next .seealso: GridGetBoundaryStart(), MeshGetBoundaryNext @*/ int GridGetBoundaryNext(Grid grid, int bd, int field, PetscTruth ghost, int *node, int *nclass) { int f; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidIntPointer(node); PetscValidIntPointer(nclass); if (grid->setupcalled == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Grid must be setup first"); /* Look for field in default class map */ for(f = 0; f < grid->cm->numFields; f++) { if (grid->cm->fields[f] == field) break; } if (f == grid->cm->numFields) { SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Field %d not present in the default class map", field); } ierr = (*grid->ops->getboundarynext)(grid, bd, field, ghost, grid->cm, node, nclass); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridRefineMesh" /*@ GridRefineMesh - Refine a grid based on area constraints. Collective on Grid Input Parameters: + grid - The initial grid - area - A function which gives an area constraint when evaluated inside an element Output Parameter: . newgrid - The refined grid Level: advanced .keywords: grid, refinement .seealso: MeshRefine(), GridReformMesh() @*/ int GridRefineMesh(Grid grid, PointFunction area, Grid *newgrid) { Mesh newMesh; int field; int ierr; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(newgrid); ierr = MeshRefine(grid->mesh, area, &newMesh); CHKERRQ(ierr); ierr = GridCreate(newMesh, newgrid); CHKERRQ(ierr); ierr = MeshDestroy(newMesh); CHKERRQ(ierr); for(field = 0; field < grid->numFields; field++) { ierr = GridAddField(*newgrid, grid->fields[field].name, grid->fields[field].discType, grid->fields[field].numComp, 0);CHKERRQ(ierr); } ierr = GridSetType(*newgrid, grid->type_name); CHKERRQ(ierr); /* Copy registered operators in discretizations */ /* Set parent: (*newgrid)->gridparent = grid; */ PetscFunctionReturn(0); } /*------------------------------------------ Distributed Data Functions ---------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridGlobalToLocalGeneral" /*@ GridGlobalToLocalGeneral - Scatters values including ghost variables from the given global vector into the given ghost vector. Collective on Grid Input Parameters: + grid - The grid . vec - The global vector . mode - The insertion mode, either INSERT_VALUES or ADD_VALUES - scatter - The scatter Output Parameters: . ghostVec - The local vector Level: intermediate .keywords: grid, scatter, global, local .seealso: GridGlobalToLocal(), GridLocalToGlobal(), GVecGlobalToLocal() @*/ int GridGlobalToLocalGeneral(Grid grid, GVec vec, Vec ghostVec, InsertMode mode, VecScatter scatter) { int ierr; PetscFunctionBegin; if (mode == ADD_VALUES) SETERRQ(PETSC_ERR_SUP, "No support for adding ghost values"); ierr = VecScatterBegin(vec, ghostVec, mode, SCATTER_FORWARD, scatter); CHKERRQ(ierr); ierr = VecScatterEnd(vec, ghostVec, mode, SCATTER_FORWARD, scatter); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGlobalToLocal" /*@ GridGlobalToLocal - Scatters values including ghost variables from the given global vector into the local grid ghost vector. Collective on Grid Input Parameters: + grid - The grid . mode - The insertion mode, either INSERT_VALUES or ADD_VALUES - vec - The grid vector Level: intermediate .keywords: grid, scatter, global, local .seealso: GridGlobalToLocalGeneral(), GridLocalToGlobal(), GVecGlobalToLocal() @*/ int GridGlobalToLocal(Grid grid, InsertMode mode, GVec vec) { Grid g; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidHeaderSpecific(vec, VEC_COOKIE); ierr = GVecGetGrid(vec, &g); CHKERRQ(ierr); if (grid != g) SETERRQ(PETSC_ERR_ARG_INCOMP, "Vector does not arise from this grid"); ierr = GridGlobalToLocalGeneral(grid, vec, grid->ghostVec, mode, grid->ghostScatter); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridLocalToGlobal" /*@ GridLocalToGlobal - Scatters local processor values including ghost variables from the local grid ghost vector into the given global vector. Collective on Grid Input Parameters: + grid - The grid . mode - The insertion mode, either SET_VALUES or ADD_VALUES - vec - The grid vector Level: intermediate .seealso: GridGlobalToLocal(), GVecLocalToGlobal() @*/ int GridLocalToGlobal(Grid grid, InsertMode mode, GVec vec) { Grid g; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidHeaderSpecific(vec, VEC_COOKIE); ierr = GVecGetGrid(vec, &g); CHKERRQ(ierr); if (grid != g) SETERRQ(PETSC_ERR_ARG_WRONG, "Vector does not arise from this grid"); ierr = VecScatterBegin(grid->ghostVec, vec, mode, SCATTER_REVERSE, grid->ghostScatter); CHKERRQ(ierr); ierr = VecScatterEnd(grid->ghostVec, vec, mode, SCATTER_REVERSE, grid->ghostScatter); CHKERRQ(ierr); PetscFunctionReturn(0); } /*---------------------------------------------- InterGrid Functions ------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridInterpolateElementVec" /*@ GridInterpolateElementVec - Interpolates a given element vector from one discretization to another. Not collective Input Parameters: + grid - The original grid . field - The original field . vec - The original element vector . newGrid - The grid defining the new element vector - newField - The new field Output Parameter: . newVec - The interpolated element vector Level: advanced .keywords: grid, interpolation .seealso: ElementVecCreate() @*/ int GridInterpolateElementVec(Grid grid, int field, ElementVec vec, Grid newGrid, int newField, ElementVec newVec) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidHeaderSpecific(newGrid, GRID_COOKIE); PetscValidHeaderSpecific(vec, ELEMENT_VEC_COOKIE); PetscValidHeaderSpecific(newVec, ELEMENT_VEC_COOKIE); GridValidField(grid, field); GridValidField(newGrid, newField); if (grid->fields[field].numComp != newGrid->fields[newField].numComp) { SETERRQ2(PETSC_ERR_ARG_INCOMP, "Fields have a different number of components %d != %d", grid->fields[field].numComp, newGrid->fields[newField].numComp); } ierr = DiscretizationInterpolateElementVec(grid->fields[field].disc, vec, newGrid->fields[newField].disc, newVec);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridCreateRestriction" /*@ GridCreateRestriction - Generates restriction matrix between two grids. Input Parameters: + vf - The fine grid - vc - The coarse grid Output Parameter: . gmat - A sparse matrix containing restriction Level: advanced @*/ int GridCreateRestriction(Grid vf, Grid vc, GMat *gmat) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(vf, GRID_COOKIE); PetscValidHeaderSpecific(vc, GRID_COOKIE); PetscValidPointer(gmat); if (vf->setupcalled == PETSC_FALSE) { ierr = GridSetUp(vf); CHKERRQ(ierr); } if (vc->setupcalled == PETSC_FALSE) { ierr = GridSetUp(vc); CHKERRQ(ierr); } if (!vf->ops->gridcreaterestriction) SETERRQ(PETSC_ERR_SUP, " "); ierr = (*vf->ops->gridcreaterestriction)(vf, vc, gmat); CHKERRQ(ierr); PetscFunctionReturn(0); } /*---------------------------------------------- Assembly Functions -------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridLocalToElementGeneral" /*@ GridLocalToElementGeneral - Scatters values including ghost variables from the given ghost vector into the given element vector. Not collective Input Parameters: + grid - The grid . ghostVec - The vector of values (including ghsot points) . reduceVec - The vector of boundary values . reduceSystem - The flag for reducing boundary conditions . reduceElement - The flag for putting boundary values in the element vector - vec - The element vector WARNING: Make sure that the indices in the element vector are local indices. Note: If reduceSystem and reduceElement are PTESC_TRUE, then boundary values are placed in vec. If reduceElement is PETSC_FALSE, then zero is used instead. Level: advanced .keywords: grid, element, scatter .seealso: GridLocalToElement(), GridGlobalToLocal() @*/ int GridLocalToElementGeneral(Grid grid, Vec ghostVec, Vec reduceVec, PetscTruth reduceSystem, PetscTruth reduceElement, ElementVec vec) { PetscScalar *array, *reduceArray; PetscScalar alpha = -grid->reduceAlpha; int elemSize = vec->reduceSize; int numOverlapVars, numReduceOverlapVars; int var; int ierr; PetscFunctionBegin; ierr = VecGetSize(ghostVec, &numOverlapVars); CHKERRQ(ierr); ierr = VecGetArray(ghostVec, &array); CHKERRQ(ierr); if (reduceSystem == PETSC_TRUE) { ierr = VecGetSize(reduceVec, &numReduceOverlapVars); CHKERRQ(ierr); ierr = VecGetArray(reduceVec, &reduceArray); CHKERRQ(ierr); for(var = 0; var < elemSize; var++) { #ifdef PETSC_USE_BOPT_g if (vec->indices[var] >= numOverlapVars) { SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Invalid index %d into ghost vector should be in [0,%d)", vec->indices[var], numOverlapVars); } #endif if (vec->indices[var] < 0) { if (reduceElement == PETSC_TRUE) { #ifdef PETSC_USE_BOPT_g if (-(vec->indices[var]+1) >= numReduceOverlapVars) { SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Invalid index %d into ghost vector should be in [0,%d)", -(vec->indices[var]+1), numReduceOverlapVars); } #endif vec->array[var] = alpha*reduceArray[-(vec->indices[var]+1)]; } else { vec->array[var] = 0.0; } } else { vec->array[var] = array[vec->indices[var]]; } } ierr = VecRestoreArray(reduceVec, &reduceArray); CHKERRQ(ierr); } else { for(var = 0; var < elemSize; var++) { #ifdef PETSC_USE_BOPT_g if ((vec->indices[var] < 0) || (vec->indices[var] >= numOverlapVars)) { SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Invalid index %d into ghost vector should be in [0,%d)", vec->indices[var], numOverlapVars); } #endif vec->array[var] = array[vec->indices[var]]; } } ierr = VecRestoreArray(ghostVec, &array); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridLocalToElement" /*@ GridLocalToElement - Scatters values including ghost variables from the local grid ghost vector into the given element vector. Not collective Input Parameters: + grid - The grid - vec - The element vector WARNING: Make sure that the indices in the element vector are local indices. Level: advanced .keywords: grid, element, scatter .seealso: GridLocalToElementGeneral(), GridGlobalToLocal() @*/ int GridLocalToElement(Grid grid, ElementVec vec) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidHeaderSpecific(vec, ELEMENT_VEC_COOKIE); ierr = GridLocalToElementGeneral(grid, grid->ghostVec, grid->bdReduceVecCur, grid->reduceSystem, grid->reduceElement, vec);CHKERRQ(ierr); PetscFunctionReturn(0); } /*-------------------------------------------- Evaluation Functions -------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridEvaluateRhs" /*@C GridEvaluateRhs - This function constructs the weak form of the functions and operators specified with GridSetRhsFunction(), GridSetRhsOperator(), and GridSetRhsNonlinearOperator(). Collective on Grid Input Parameters: + grid - The grid . x - The current iterate - ctx - The optional user context Output Parameters . f - The function value Level: beginner .keywords: grid, rhs .seealso: GridEvaluateSystemMatrix(), GridSetRhsFunction(), GridSetRhsOperator(), GridSetRhsNonlinearOperator() @*/ int GridEvaluateRhs(Grid grid, GVec x, GVec f, PetscObject ctx) { Mesh mesh; MeshMover mover; Grid meshVelGrid; PetscScalar alpha = -grid->reduceAlpha; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); /* x can be PETSC_NULL */ PetscValidHeaderSpecific(f, VEC_COOKIE); /* Setup mesh velocity calculation for ALE variables */ if (grid->ALEActive == PETSC_TRUE) { ierr = GridGetMesh(grid, &mesh); CHKERRQ(ierr); ierr = MeshGetMover(mesh, &mover); CHKERRQ(ierr); ierr = MeshMoverGetVelocityGrid(mover, &meshVelGrid); CHKERRQ(ierr); ierr = GridSetUp(meshVelGrid); CHKERRQ(ierr); } /* Calculate Rhs */ ierr = (*grid->ops->gridevaluaterhs)(grid, x, f, ctx); CHKERRQ(ierr); /* Add boundary values */ if ((grid->reduceSystem == PETSC_TRUE) && (grid->reduceElement == PETSC_FALSE) && (grid->bdReduceMat != PETSC_NULL)) { /* Calculate contribution of constrained variables to the Rhs */ ierr = MatMult(grid->bdReduceMat, grid->bdReduceVecCur, grid->reduceVec); CHKERRQ(ierr); /* Update Rhs */ ierr = VecAXPY(&alpha, grid->reduceVec, f); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridEvaluateRhsFunction" /*@C GridEvaluateRhsFunction This function constructs the weak form of the functions specified with GridSetRhsFunction(). Collective on Grid Input Parameters: + grid - The grid . x - The current iterate - ctx - The optional user context Output Parameter: . f - The function value Level: beginner .keywords grid, rhs, function .seealso GridEvaluateRhsOperator(), GridEvaluateRhs(), GridEvaluateSystemMatrix() @*/ int GridEvaluateRhsFunction(Grid grid, GVec x, GVec f, void *ctx) { PetscTruth oldLinear = grid->activeOpTypes[1]; PetscTruth oldNonlinear = grid->activeOpTypes[2]; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); /* x can be PETSC_NULL */ PetscValidHeaderSpecific(f, VEC_COOKIE); /* Turn off computation of operators */ grid->activeOpTypes[1] = PETSC_FALSE; grid->activeOpTypes[2] = PETSC_FALSE; ierr = GridEvaluateRhs(grid, x, f, (PetscObject) ctx); CHKERRQ(ierr); grid->activeOpTypes[1] = oldLinear; grid->activeOpTypes[2] = oldNonlinear; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridEvaluateRhsOperator" /*@C GridEvaluateRhsOperator This function constructs the weak form of the operators specified with GridSetRhsOperator(). Collective on Grid Input Parameters: + grid - The grid . x - The current iterate . linear - The flag for including linear operators . nonlinear - The flag for including nonlinear operators - ctx - The optional user context Output Parameter: . f - The operator value Level: beginner .keywords grid, rhs, operator .seealso GridEvaluateRhsFunction(), GridEvaluateRhs(), GridEvaluateSystemMatrix() @*/ int GridEvaluateRhsOperator(Grid grid, GVec x, GVec f, PetscTruth linear, PetscTruth nonlinear, void *ctx) { PetscTruth oldFunction = grid->activeOpTypes[0]; PetscTruth oldLinear = grid->activeOpTypes[1]; PetscTruth oldNonlinear = grid->activeOpTypes[2]; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); /* x can be PETSC_NULL */ PetscValidHeaderSpecific(f, VEC_COOKIE); /* Turn off computation of operators */ grid->activeOpTypes[0] = PETSC_FALSE; grid->activeOpTypes[1] = linear; grid->activeOpTypes[2] = nonlinear; ierr = GridEvaluateRhs(grid, x, f, (PetscObject) ctx); CHKERRQ(ierr); grid->activeOpTypes[0] = oldFunction; grid->activeOpTypes[1] = oldLinear; grid->activeOpTypes[2] = oldNonlinear; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridEvaluateSystemMatrix" /*@C GridEvaluateSystemMatrix This function constructs the weak form of the linear operators specified with GridSetMatOperator(). Collective on Grid Input Parameters: + grid - The grid . x - The current iterate - ctx - The optional user context Output Parameter: . f - The function value Level: beginner .keywords: grid, matrix .seealso: GridEvaluateRhs() @*/ int GridEvaluateSystemMatrix(Grid grid, GVec x, GMat *J, GMat *M, MatStructure *flag, PetscObject ctx) { Mesh mesh; MeshMover mover; Grid meshVelGrid; #ifdef PETSC_HAVE_MATHEMATICA PetscTruth opt; #endif int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); if (grid->ALEActive == PETSC_TRUE) { ierr = GridGetMesh(grid, &mesh); CHKERRQ(ierr); ierr = MeshGetMover(mesh, &mover); CHKERRQ(ierr); ierr = MeshMoverGetVelocityGrid(mover, &meshVelGrid); CHKERRQ(ierr); ierr = GridSetUp(meshVelGrid); CHKERRQ(ierr); } if (grid->isMatrixFree == PETSC_TRUE) { /* Should probably set this as the matrix context, but there is no MatShellSetContext() */ grid->matrixFreeArg = x; grid->matrixFreeContext = ctx; } else { ierr = (*grid->ops->gridevaluatesystemmatrix)(grid, x, J, M, flag, ctx); CHKERRQ(ierr); } #ifdef PETSC_HAVE_MATHEMATICA /* Debugging */ ierr = PetscOptionsHasName(PETSC_NULL, "-trace_grid_math", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = ViewerMathematicaSetName(PETSC_VIEWER_MATHEMATICA_WORLD, "jac"); CHKERRQ(ierr); ierr = MatView(*J, PETSC_VIEWER_MATHEMATICA_WORLD); CHKERRQ(ierr); ierr = ViewerMathematicaClearName(PETSC_VIEWER_MATHEMATICA_WORLD); CHKERRQ(ierr); } #endif PetscFunctionReturn(0); } #if 0 /*----------------------------------------- Parallel Communication Functions ----------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridGhostExchange" /*@C GridGhostExchange - This functions transfers data between local and ghost storage without a predefined mapping. Collective on MPI_Comm Input Parameters: + comm - The communicator . numGhosts - The number of ghosts in this domain . ghostProcs - The processor from which to obtain each ghost . ghostIndices - The global index for each ghost . dataType - The type of the variables . firstVar - The first variable on each processor . addv - The insert mode, INSERT_VALUES or ADD_VALUES - mode - The direction of the transfer, SCATTER_FORWARD or SCATTER_REVERSE Output Paramters: + locVars - The local variable array - ghostVars - The ghost variables Note: The data in ghostVars is assumed contiguous and implicitly indexed by the order of ghostProcs and ghostIndices. The SCATTER_FORWARD mode will take the requested data from locVars and copy it to ghostVars in the order specified by ghostIndices. The SCATTER_REVERSE mode will take data from ghostVars and copy it to locVars. Level: advanced .keywords: ghost, exchange, grid .seealso: GridGlobalToLocal(), GridLocalToGlobal() @*/ int GridGhostExchange(MPI_Comm comm, int numGhosts, int *ghostProcs, int *ghostIndices, PetscDataType dataType, int *firstVar, InsertMode addv, ScatterMode mode, void *locVars, void *ghostVars) { int *numSendGhosts; /* The number of ghosts from each domain */ int *numRecvGhosts; /* The number of local variables which are ghosts in each domain */ int *sumSendGhosts; /* The prefix sums of numSendGhosts */ int *sumRecvGhosts; /* The prefix sums of numRecvGhosts */ int *offsets; /* The offset into the send array for each domain */ int totSendGhosts; /* The number of ghosts to request variables for */ int totRecvGhosts; /* The number of nodes to provide class info about */ int *sendIndices; /* The canonical indices of ghosts in this domain */ int *recvIndices; /* The canonical indices of ghosts to return variables for */ char *tempVars; /* The variables of the requested or submitted ghosts */ char *locBytes = (char *) locVars; MPI_Datatype MPIType; int typeSize; #ifdef PETSC_USE_BOPT_g int numLocVars; #endif int numProcs, rank; int proc, ghost, locIndex, byte; int ierr; PetscFunctionBegin; /* Initialize communication */ ierr = MPI_Comm_size(comm, &numProcs); CHKERRQ(ierr); ierr = MPI_Comm_rank(comm, &rank); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &numSendGhosts); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &numRecvGhosts); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &sumSendGhosts); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &sumRecvGhosts); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &offsets); CHKERRQ(ierr); ierr = PetscMemzero(numSendGhosts, numProcs * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemzero(numRecvGhosts, numProcs * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemzero(sumSendGhosts, numProcs * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemzero(sumRecvGhosts, numProcs * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemzero(offsets, numProcs * sizeof(int)); CHKERRQ(ierr); #ifdef PETSC_USE_BOPT_g numLocVars = firstVar[rank+1] - firstVar[rank]; #endif /* Get number of ghosts needed from each processor */ for(ghost = 0; ghost < numGhosts; ghost++) numSendGhosts[ghostProcs[ghost]]++; /* Get number of ghosts to provide variables for */ ierr = MPI_Alltoall(numSendGhosts, 1, MPI_INT, numRecvGhosts, 1, MPI_INT, comm); CHKERRQ(ierr); for(proc = 1; proc < numProcs; proc++) { sumSendGhosts[proc] = sumSendGhosts[proc-1] + numSendGhosts[proc-1]; sumRecvGhosts[proc] = sumRecvGhosts[proc-1] + numRecvGhosts[proc-1]; offsets[proc] = sumSendGhosts[proc]; } totSendGhosts = sumSendGhosts[numProcs-1] + numSendGhosts[numProcs-1]; totRecvGhosts = sumRecvGhosts[numProcs-1] + numRecvGhosts[numProcs-1]; if (numGhosts != totSendGhosts) { SETERRQ2(PETSC_ERR_PLIB, "Invalid number of ghosts %d in send, should be %d", totSendGhosts, numGhosts); } ierr = PetscDataTypeGetSize(dataType, &typeSize); CHKERRQ(ierr); if (totSendGhosts) { sendIndices = (int *) PetscMalloc(totSendGhosts * sizeof(int)); CHKERRQ(sendIndices); } if (totRecvGhosts) { recvIndices = (int *) PetscMalloc(totRecvGhosts * sizeof(int)); CHKERRQ(recvIndices); tempVars = (char *) PetscMalloc(totRecvGhosts * typeSize); CHKERRQ(tempVars); } /* Must order ghosts by processor */ for(ghost = 0; ghost < numGhosts; ghost++) sendIndices[offsets[ghostProcs[ghost]]++] = ghostIndices[ghost]; /* Get canonical indices of ghosts to provide variables for */ ierr = MPI_Alltoallv(sendIndices, numSendGhosts, sumSendGhosts, MPI_INT, recvIndices, numRecvGhosts, sumRecvGhosts, MPI_INT, comm); CHKERRQ(ierr); switch(mode) { case SCATTER_FORWARD: /* Get ghost variables */ if (addv == INSERT_VALUES) { for(ghost = 0; ghost < totRecvGhosts; ghost++) { locIndex = recvIndices[ghost] - firstVar[rank]; #ifdef PETSC_USE_BOPT_g if ((locIndex < 0) || (locIndex >= numLocVars)) { SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars); } #endif for(byte = 0; byte < typeSize; byte++) tempVars[ghost*typeSize+byte] = locBytes[locIndex*typeSize+byte]; } } else { for(ghost = 0; ghost < totRecvGhosts; ghost++) { locIndex = recvIndices[ghost] - firstVar[rank]; #ifdef PETSC_USE_BOPT_g if ((locIndex < 0) || (locIndex >= numLocVars)) { SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars); } #endif for(byte = 0; byte < typeSize; byte++) tempVars[ghost*typeSize+byte] += locBytes[locIndex*typeSize+byte]; } } /* Communicate local variables to ghost storage */ ierr = PetscDataTypeToMPIDataType(dataType, &MPIType); CHKERRQ(ierr); ierr = MPI_Alltoallv(tempVars, numRecvGhosts, sumRecvGhosts, MPIType, ghostVars, numSendGhosts, sumSendGhosts, MPIType, comm); CHKERRQ(ierr); break; case SCATTER_REVERSE: /* Communicate ghost variables to local storage */ ierr = PetscDataTypeToMPIDataType(dataType, &MPIType); CHKERRQ(ierr); ierr = MPI_Alltoallv(ghostVars, numSendGhosts, sumSendGhosts, MPIType, tempVars, numRecvGhosts, sumRecvGhosts, MPIType, comm); CHKERRQ(ierr); /* Get ghost variables */ if (addv == INSERT_VALUES) { for(ghost = 0; ghost < totRecvGhosts; ghost++) { locIndex = recvIndices[ghost] - firstVar[rank]; #ifdef PETSC_USE_BOPT_g if ((locIndex < 0) || (locIndex >= numLocVars)) { SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars); } #endif for(byte = 0; byte < typeSize; byte++) locBytes[locIndex*typeSize+byte] = tempVars[ghost*typeSize+byte]; } } else { /* There must be a better way to do this -- Ask Bill */ if (typeSize == sizeof(int)) { int *tempInt = (int *) tempVars; int *locInt = (int *) locVars; for(ghost = 0; ghost < totRecvGhosts; ghost++) { locIndex = recvIndices[ghost] - firstVar[rank]; #ifdef PETSC_USE_BOPT_g if ((locIndex < 0) || (locIndex >= numLocVars)) { SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars); } #endif locInt[locIndex] += tempInt[ghost]; } } else if (typeSize == sizeof(long int)) { long int *tempLongInt = (long int *) tempVars; long int *locLongInt = (long int *) locVars; for(ghost = 0; ghost < totRecvGhosts; ghost++) { locIndex = recvIndices[ghost] - firstVar[rank]; #ifdef PETSC_USE_BOPT_g if ((locIndex < 0) || (locIndex >= numLocVars)) { SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars); } #endif locLongInt[locIndex] += tempLongInt[ghost]; } } else { for(ghost = 0; ghost < totRecvGhosts; ghost++) { locIndex = recvIndices[ghost] - firstVar[rank]; #ifdef PETSC_USE_BOPT_g if ((locIndex < 0) || (locIndex >= numLocVars)) { SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid ghost index %d, not in [0,%d)", locIndex, numLocVars); } #endif for(byte = 0; byte < typeSize; byte++) locBytes[locIndex*typeSize+byte] += tempVars[ghost*typeSize+byte]; } } } break; default: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid scatter mode %d", mode); } /* Cleanup */ ierr = PetscFree(numSendGhosts); CHKERRQ(ierr); ierr = PetscFree(numRecvGhosts); CHKERRQ(ierr); ierr = PetscFree(sumSendGhosts); CHKERRQ(ierr); ierr = PetscFree(sumRecvGhosts); CHKERRQ(ierr); ierr = PetscFree(offsets); CHKERRQ(ierr); if (totSendGhosts) { ierr = PetscFree(sendIndices); CHKERRQ(ierr); } if (totRecvGhosts) { ierr = PetscFree(recvIndices); CHKERRQ(ierr); ierr = PetscFree(tempVars); CHKERRQ(ierr); } PetscFunctionReturn(0); } #endif