#ifdef PETSC_RCS_HEADER static char vcid[] = "$Id: gvec.c,v 1.48 2000/07/16 23:20:04 knepley Exp $"; #endif /* This file provides routines for grid vectors (vectors that are associated with grids, possibly with multiple degrees of freedom per node). */ #include "petscsles.h" /* For ALE Mesh Support */ #include "src/gvec/gvecimpl.h" /*I "gvec.h" I*/ /* Logging support */ int GVEC_EvaluateFunction, GVEC_EvaluateFunctionCollective, GVEC_EvaluateRhs, GVEC_EvaluateJacobian; int GVEC_SetBoundary, GVEC_InterpolateField, GVEC_InterpolateFieldBatch, GVEC_InterpolateFieldBatchParallel; int GVEC_InterpolateFieldBatchCalc; extern int VecCreate_MPI_Private(Vec, int, const PetscScalar [], PetscMap); /*--------------------------------------------- General Vector Functions ---------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GVecDestroy" /*@C GVecDestroy - Destroys a grid vector. Input Parameters: . v - the vector Level: beginner .keywords: grid vector, destroy .seealso: VecDestroy(), GVecDuplicate() @*/ int GVecDestroy(GVec v) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); if (--v->refct > 0) PetscFunctionReturn(0); ierr = VecDestroy(v); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecViewFromOptions" /*@ GVecViewFromOptions - This function visualizes the vector based upon user options. Collective on GVec Input Parameter: . gvec - The vector Level: intermediate .keywords: GVec, view, options, database .seealso: GVecSetFromOptions(), GVecView() @*/ int GVecViewFromOptions(GVec gvec, 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(gvec->prefix, "-gvec_view", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = PetscOptionsGetString(gvec->prefix, "-gvec_view", typeName, 1024, &opt); CHKERRQ(ierr); ierr = PetscStrlen(typeName, &len); CHKERRQ(ierr); if (len > 0) { ierr = PetscViewerCreate(gvec->comm, &viewer); CHKERRQ(ierr); ierr = PetscViewerSetType(viewer, typeName); CHKERRQ(ierr); ierr = PetscOptionsGetString(gvec->prefix, "-gvec_view_file", fileName, 1024, &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = PetscViewerSetFilename(viewer, fileName); CHKERRQ(ierr); } else { ierr = PetscViewerSetFilename(viewer, gvec->name); CHKERRQ(ierr); } ierr = GVecView(gvec, viewer); CHKERRQ(ierr); ierr = PetscViewerFlush(viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(viewer); CHKERRQ(ierr); } else { ierr = GVecView(gvec, PETSC_NULL); CHKERRQ(ierr); } } ierr = PetscOptionsHasName(gvec->prefix, "-gvec_view_draw", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = PetscViewerDrawOpen(gvec->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) gvec); CHKERRQ(ierr) ; titleStr = gvec->name; } ierr = PetscDrawSetTitle(draw, titleStr); CHKERRQ(ierr); ierr = GVecView(gvec, viewer); CHKERRQ(ierr); ierr = PetscViewerFlush(viewer); CHKERRQ(ierr); ierr = PetscDrawPause(draw); CHKERRQ(ierr); ierr = PetscViewerDestroy(viewer); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecView" /*@ GVecView - Views a grid vector. Input Parameters: . v - The grid vector . viewer - [Optional] A visualization context Notes: GVecView() supports the same viewers as VecView(). The only difference is that for all multiprocessor cases, the output vector employs the natural ordering of the grid, so it many cases this corresponds to the ordering that would have been used for the uniprocessor case. The available visualization contexts include $ PETSC_VIEWER_STDOUT_SELF - standard output (default) $ PETSC_VIEWER_STDOUT_WORLD - synchronized standard $ output where only the first processor opens $ the file. All other processors send their $ data to the first processor to print. The user can open alternative vistualization contexts with $ PetscViewerFileOpenASCII() - output vector to a specified file $ PetscViewerFileOpenBinary() - output in binary to a $ specified file; corresponding input uses VecLoad() $ PetscViewerDrawOpenX() - output vector to an X window display $ DrawLGCreate() - output vector as a line graph to an X window display $ PetscViewerMatlabOpen() - output vector to Matlab viewer $ PetscViewerMathematicaOpen() - output vector to Mathematica viewer Level: beginner .keywords: view, visualize, output, print, write, draw .seealso: VecView() @*/ int GVecView(GVec v, PetscViewer viewer) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); if (!viewer) { viewer = PETSC_VIEWER_STDOUT_SELF; } else { PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE); } ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = (*grid->ops->gvecview)(v, viewer); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecSerialize" /*@ GVecSerialize - This function stores or recreates a grid vector using a viewer for a binary file. Input Parameters: . viewer - The viewer context . store - This flag is PETSC_TRUE is data is being written, otherwise it will be read Output Parameter: . v - The grid vector Level: beginner .keywords: grid vector, serialize .seealso: GridSerialize() @*/ int GVecSerialize(Grid grid, GVec *v, PetscViewer viewer, PetscTruth store) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE); PetscValidPointer(v); ierr = VecSerialize(grid->comm, v, viewer, store); CHKERRQ(ierr); /* Newly created object should have grid as a parent */ if (store == PETSC_FALSE) { ierr = PetscObjectCompose((PetscObject) *v, "Grid", (PetscObject) grid); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecDuplicate" /*@C GVecDuplicate - Duplicates a grid vector. Input Parameters: . v - the vector Level: beginner .keywords: grid vector, destroy .seealso: VecDuplicate(), GVecDestroy() @*/ int GVecDuplicate(GVec v, GVec *newvec) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); ierr = VecDuplicate(v, newvec); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecGetGrid" /*@ GVecGetGrid - This function returns the grid from a grid vector. Not collective Input Parameter: . v - The vector Output Parameter: . grid - The grid Level: intermediate .keywords: grid vector, get, grid .seealso: GVecGetOrder(), GridGetMesh(), GMatGetGrid() @*/ int GVecGetGrid(GVec v, Grid *grid) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidPointer(grid); ierr = PetscObjectQuery((PetscObject) v, "Grid", (PetscObject *) grid); CHKERRQ(ierr); PetscValidHeaderSpecific(*grid, GRID_COOKIE); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecGetOrder" /*@ GVecGetOrder - This function returns the variable ordering from a grid vector. Not collective Input Parameter: . v - The vector Output Parameter: . order - The variable ordering Level: intermediate .keywords: grid vector, get, variable ordering .seealso: GVecGetGrid(), GridGetMesh(), GMatGetGrid() @*/ int GVecGetOrder(GVec v, VarOrdering *order) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidPointer(order); ierr = PetscObjectQuery((PetscObject) v, "Order", (PetscObject *) order); CHKERRQ(ierr); PetscValidHeaderSpecific(*order, VAR_ORDER_COOKIE); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecScatterCreate" /*@C GVecScatterCreate - This function creates a scatter from a vector to an embedded vector. Collective on GVec Input Parameters: + vec - A vector - embedVec - A vector constructed from an embedded ordering Output Parameter: . scatter - The scatter Level: beginner .keywords variable ordering, vector scatter, embedding .seealso VecScatterCreate() @*/ int GVecScatterCreate(GVec vec, GVec embedVec, VecScatter *scatter) { Grid grid; VarOrdering order, embedOrder; FieldClassMap map, embedMap; int f, field, nclass; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(vec, VEC_COOKIE); PetscValidHeaderSpecific(embedVec, VEC_COOKIE); PetscValidPointer(scatter); ierr = GVecGetGrid(vec, &grid); CHKERRQ(ierr); ierr = GVecGetOrder(vec, &order); CHKERRQ(ierr); ierr = GVecGetOrder(embedVec, &embedOrder); CHKERRQ(ierr); ierr = VarOrderingGetClassMap(order, &map); CHKERRQ(ierr); ierr = VarOrderingGetClassMap(embedOrder, &embedMap); CHKERRQ(ierr); if (map->numFields < embedMap->numFields) { SETERRQ2(PETSC_ERR_ARG_INCOMP, "The embedded ordering cannot have fewer fields than the original (%d < %d)", map->numFields, embedMap->numFields); } if (map->numClasses != embedMap->numClasses) { SETERRQ2(PETSC_ERR_ARG_INCOMP, "The embedded ordering must have the same number of classes as the original (%d < %d)", map->numClasses, embedMap->numClasses); } for(nclass = 0; nclass < map->numClasses; nclass++) { if (embedMap->classSizes[nclass] > 0) { if (map->classSizes[nclass] <= 0) { SETERRQ1(PETSC_ERR_ARG_INCOMP, "Mapping not an embedding: class %d is larger than the original", nclass); } } for(f = 0; f < embedMap->numFields; f++) { field = map->fields[f]; if ((order->localStart[field][nclass] < 0) && (embedOrder->localStart[field][nclass] >= 0)) { SETERRQ2(PETSC_ERR_ARG_INCOMP, "Mapping not an embedding: field %d not present in class %d in the original", field, nclass); } } } ierr = GridSetUp(grid); CHKERRQ(ierr); ierr = (*grid->ops->gridcreatevarscatter)(grid, order, vec, embedOrder, embedVec, scatter); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecGetLocalWorkGVec" /*@ GVecGetLocalWorkGVec - Creates a new GVec to contain the local + ghost values portion of a GVec. Input Parameter: . v - the grid vector Output Parameters: . lv - the local vector Level: intermediate .seealso: GVecLocalToGlobal(), GVecGlobalToLocal(), GVecRestoreLocalWorkGVec(), GVecGetWorkGVec(), GVecRestoreWorkGVec() @*/ int GVecGetLocalWorkGVec(GVec v, GVec *lv) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidPointer(lv); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = (*grid->ops->gvecgetlocalworkgvec)(v, lv); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecRestoreLocalWorkGVec" /*@ GVecRestoreLocalWorkGVec - Input Parameter: . v - the grid vector Output Parameters: . lv - the local vector Level: intermediate .seealso: GVecLocalToGlobal(), GVecGlobalToLocal(), GVecGetLocalWorkGVec(), GVecGetWorkGVec(), GVecRestoreWorkGVec() @*/ int GVecRestoreLocalWorkGVec(GVec v, GVec *lv) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidPointer(lv); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = (*grid->ops->gvecrestorelocalworkgvec)(v, lv); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecGetWorkGVec" /*@ GVecGetWorkGVec - Get an identical work grid vector. Input Parameter: . v - The grid vector Output Parameters: . wv - The work vector Level: intermediate .seealso: GVecLocalToGlobal(), GVecGlobalToLocal(), GVecRestoreWorkGVec(), GVecGetLocalWorkGVec(), GVecRestoreLocalWorkGVec() @*/ int GVecGetWorkGVec(GVec v, GVec *wv) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidPointer(wv); ierr = VecDuplicate(v, wv); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecRestoreWorkGVec" /*@ GVecRestoreWorkGVec - Restores the work vector. Input Parameter: . v - The grid vector Output Parameters: . wv - The work vector Level: intermediate .seealso: GVecLocalToGlobal(), GVecGlobalToLocal(), GVecGetWorkGVec(), GVecGetLocalWorkGVec(), GVecRestoreLocalWorkGVec() @*/ int GVecRestoreWorkGVec(GVec v, GVec *wv) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidPointer(wv); ierr = VecDestroy(*wv); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecGlobalToLocal" /*@ GVecGlobalToLocal - Copies a global vector including ghost points to a local work vector. Input Parameter: . v - the grid vector . mode - either SET_VALUES or ADD_VALUES Output Parameters: . lv - the local vector Level: intermediate .seealso: GVecLocalToGlobal(), GVecCreateLocalGVec() @*/ int GVecGlobalToLocal(GVec v, InsertMode mode, GVec lv) { int ierr; Grid grid; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidHeaderSpecific(lv, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = (*grid->ops->gvecglobaltolocal)(v, mode, lv); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecLocalToGlobal" /*@ GVecLocalToGlobal - Copies a local vector including ghost points to its global representation. Input Parameter: . v - the local grid vector . mode - either SET_VALUES or ADD_VALUES Output Parameters: . lv - the global vector Level: intermediate .seealso: GVecGlobalToLocal(), GVecCreateLocalGVec() @*/ int GVecLocalToGlobal(GVec v,InsertMode mode,GVec lv) { int ierr; Grid grid; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidHeaderSpecific(lv, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = (*grid->ops->gveclocaltoglobal)(v, mode, lv); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateFunction" /*@ GVecEvaluateFunction - Evaluates a function over the specified fields on the locations defined by the underlying grid and its discretization. Input Parameters: + v - The grid vector . numFields - The number of fields to evaluate . fields - The fields . f - The user provided function . alpha - The scalar multiplier - ctx - [Optional] A user provided context for the function The function f should return ordered first by node and then by component. For instance, if superscripts denote components and subscripts denote nodes: v^0_0 v^1_0 v^2_0 v^0_1 v^1_1 v^2_1 ... v^0_n v^1_n v^2_n This is the standard for PointFunctions. Level: intermediate .seealso: GVecEvaluateFunctionGalerkin,GMatEvaluateOperatorGalerkin @*/ int GVecEvaluateFunction(GVec v, int numFields, int *fields, PointFunction f, PetscScalar alpha, void *ctx) { Grid grid; VarOrdering order; int size; int fieldIdx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = PetscLogEventBegin(GVEC_EvaluateFunction, v, grid, 0, 0); CHKERRQ(ierr); ierr = VecGetSize(v, &size); CHKERRQ(ierr); for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) if ((fields[fieldIdx] < 0) || (fields[fieldIdx] > grid->numFields)) { SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid field number %d", fields[fieldIdx]); } /* Should check that fields are in the local order */ if ((grid->isConstrained) && (size == grid->constraintOrder->numVars)) { ierr = VarOrderingCreateSubset(grid->constraintOrder, numFields, fields, PETSC_FALSE, &order); CHKERRQ(ierr); /* This is a kludge to get around the size check since the new constrained vars are not formed */ order->numLocVars -= order->numLocNewVars; ierr = (*grid->ops->gvecevaluatefunction)(grid, v, order, f, alpha, ctx); CHKERRQ(ierr); order->numLocVars += order->numLocNewVars; } else if (size == grid->order->numVars) { ierr = VarOrderingCreateSubset(grid->order, numFields, fields, PETSC_FALSE, &order); CHKERRQ(ierr); ierr = (*grid->ops->gvecevaluatefunction)(grid, v, order, f, alpha, ctx); CHKERRQ(ierr); } else { SETERRQ(PETSC_ERR_ARG_INCOMP, "Vector size does not conform to grid"); } ierr = VarOrderingDestroy(order); CHKERRQ(ierr); ierr = PetscLogEventEnd(GVEC_EvaluateFunction, v, grid, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateFunctionCollective" /*@ GVecEvaluateFunctionCollective - Evaluates a function over the specified fields on the locations defined by the underlying grid and its discretization. The only difference between it and GVecEvaluateFunction(), is that each processor is guaranteed to call f at each iteration, possibly with null arguments. Input Parameters: + v - The grid vector . numFields - The number of fields to evaluate . fields - The fields . f - The user provided function . alpha - The scalar multiplier - ctx - [Optional] The user provided context for the function The function f should return ordered first by node and then by component. For instance, if superscripts denote components and subscripts denote nodes: v^0_0 v^1_0 v^2_0 v^0_1 v^1_1 v^2_1 ... v^0_n v^1_n v^2_n This is the standard for PointFunctions. Also, this function is typically used with an f that has a barrier. In this way you can guarantee that any collective operation will complete inside f. Level: intermediate .seealso: GVecEvaluateFunctionGalerkin,GMatEvaluateOperatorGalerkin @*/ int GVecEvaluateFunctionCollective(GVec v, int numFields, int *fields, PointFunction f, PetscScalar alpha, void *ctx) { Grid grid; VarOrdering order; int size; int fieldIdx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = PetscLogEventBegin(GVEC_EvaluateFunctionCollective, v, grid, 0, 0); CHKERRQ(ierr); ierr = VecGetSize(v, &size); CHKERRQ(ierr); for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) if ((fields[fieldIdx] < 0) || (fields[fieldIdx] > grid->numFields)) { SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid field number %d", fields[fieldIdx]); } /* Should check that fields are in the local order */ if ((grid->isConstrained) && (size == grid->constraintOrder->numVars)) { ierr = VarOrderingCreateSubset(grid->constraintOrder, numFields, fields, PETSC_FALSE, &order); CHKERRQ(ierr); /* This is a kludge to get around the size check since the new constrained vars are not formed */ order->numLocVars -= order->numLocNewVars; ierr = (*grid->ops->gvecevaluatefunctioncollective)(grid, v, order, f, alpha, ctx); CHKERRQ(ierr); order->numLocVars += order->numLocNewVars; } else if (size == grid->order->numVars) { ierr = VarOrderingCreateSubset(grid->order, numFields, fields, PETSC_FALSE, &order); CHKERRQ(ierr); ierr = (*grid->ops->gvecevaluatefunctioncollective)(grid, v, order, f, alpha, ctx); CHKERRQ(ierr); } else { SETERRQ(PETSC_ERR_ARG_INCOMP, "Vector size does not conform to grid"); } ierr = VarOrderingDestroy(order); CHKERRQ(ierr); ierr = PetscLogEventEnd(GVEC_EvaluateFunctionCollective, v, grid, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateFunctionRectangular" /*@ GVecEvaluateFunctionRectangular - Evaluates a function over the specified fields on the locations defined by the underlying grid and its discretization, and takes a user-defined ordering. Input Parameters: + v - The grid vector . order - The variable ordering . f - The user provided function . alpha - A scalar multiplier - ctx - An optional user provided context for the function Level: intermediate .seealso: GVecEvaluateFunctionGalerkin,GMatEvaluateOperatorGalerkin @*/ int GVecEvaluateFunctionRectangular(GVec v, VarOrdering order, PointFunction f, PetscScalar alpha, void *ctx) { Grid grid; int size; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = VecGetSize(v, &size); CHKERRQ(ierr); if (size != order->numVars) SETERRQ(PETSC_ERR_ARG_INCOMP, "Vector size does not conform to ordering"); /* Should check that fields are in the local order */ ierr = (*grid->ops->gvecevaluatefunction)(grid, v, order, f, alpha, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateFunctionGalerkin" /*@ GVecEvaluateFunctionGalerkin - Evaluates the weak form of a function over a field on the locations defined by the underlying grid and its discretization. Input Parameter: . v - The grid vector . numFields - The number of fields to evaluate . fields - The fields . f - The user provided function . alpha - A scalar multiplier . ctx - An optional user provided context for the function Level: intermediate .seealso: GVecEvaluateFunction,GMatEvaluateFunctionGalerkin @*/ int GVecEvaluateFunctionGalerkin(GVec v, int numFields, int *fields, PointFunction f, PetscScalar alpha, void *ctx) { Grid grid; int fieldIdx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidIntPointer(fields); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) { if ((fields[fieldIdx] < 0) || (fields[fieldIdx] > grid->numFields)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid field number"); } } /* Should check that fields are in the local order */ ierr = (*grid->ops->gvecevaluatefunctiongalerkin)(grid, v, numFields, fields, grid->locOrder, f, alpha, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateFunctionGalerkinCollective" /*@ GVecEvaluateFunctionGalerkinCollective - Evaluates the weak form of a function over a field on the locations defined by the underlying grid and its discretization. The only difference between it and GVecEvaluateFunctionGalerkin(), is that each processor is guaranteed to call f at each iteration, possibly with null arguments. Input Parameter: . v - The grid vector . numFields - The number of fields to evaluate . fields - The fields . f - The user provided function . alpha - A scalar multiplier . ctx - An optional user provided context for the function Level: intermediate .seealso: GVecEvaluateFunction,GMatEvaluateFunctionGalerkin @*/ int GVecEvaluateFunctionGalerkinCollective(GVec v, int numFields, int *fields, PointFunction f, PetscScalar alpha, void *ctx) { Grid grid; int fieldIdx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidIntPointer(fields); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) { if ((fields[fieldIdx] < 0) || (fields[fieldIdx] > grid->numFields)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid field number"); } } /* Should check that fields are in the local order */ ierr = (*grid->ops->gvecevaluatefunctiongalerkincollective)(grid, v, numFields, fields, grid->locOrder, f, alpha, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateBoundaryFunctionGalerkin" /*@ GVecEvaluateBoundaryFunctionGalerkin - Evaluates the weak form of a function over a field on the locations defined by the boundary of the underlying mesh and its discretization. Input Parameter: . v - The grid vector . numFields - The number of fields to evaluate . fields - The fields . f - The user provided function . alpha - A scalar multiplier . ctx - An optional user provided context for the function Level: intermediate .seealso: GVecEvaluateFunction,GMatEvaluateFunctionGalerkin @*/ int GVecEvaluateBoundaryFunctionGalerkin(GVec v, int numFields, int *fields, PointFunction f, PetscScalar alpha, void *ctx) { Grid grid; int fieldIdx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidIntPointer(fields); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) { if ((fields[fieldIdx] < 0) || (fields[fieldIdx] > grid->numFields)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid field number"); } } /* Should check that fields are in the local order */ ierr = (*grid->ops->gvecevaluateboundaryfunctiongalerkin)(grid, v, numFields, fields, grid->bdLocOrder, f, alpha, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateBoundaryFunctionGalerkinCollective" /*@ GVecEvaluateBoundaryFunctionGalerkinCollective - Evaluates the weak form of a function over a field on the locations defined by the boundary of the underlying mesh and its discretization. The only difference between it and GVecEvaluateBoundaryFunctionGalerkin(), is that each processor is guaranteed to call f at each iteration, possibly with null arguments. Input Parameter: . v - The grid vector . numFields - The number of fields to evaluate . fields - The fields . f - The user provided function . alpha - A scalar multiplier . ctx - An optional user provided context for the function Level: intermediate .seealso: GVecEvaluateFunction,GMatEvaluateFunctionGalerkin @*/ int GVecEvaluateBoundaryFunctionGalerkinCollective(GVec v, int numFields, int *fields, PointFunction f, PetscScalar alpha, void *ctx) { Grid grid; int fieldIdx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidIntPointer(fields); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) { if ((fields[fieldIdx] < 0) || (fields[fieldIdx] > grid->numFields)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid field number"); } } /* Should check that fields are in the local order */ ierr = (*grid->ops->gvecevaluateboundaryfunctiongalerkincollective)(grid, v, numFields, fields, grid->bdLocOrder, f, alpha, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateNonlinearOperatorGalerkin" /*@ GVecEvaluateNonlinearOperatorGalerkin - Evaluates the weak form of a nonlinear operator over a field on the locations defined by the underlying grid and its discretization. Collective on GVec Input Parameter: . v - The grid vector . n - The number of input vectors . vecs - The input vectors . numFields - The number of fields to evaluate . fields - The fields . op - The nonlinear operator . alpha - The scalar multiplier . isALE - The flag for ALE operators . ctx - [Optional] A user provided context for the function Level: intermediate .keywords: grid vector, nonlinear, operator, galerkin .seealso: GVecEvaluateFunctionGalerkin() @*/ int GVecEvaluateNonlinearOperatorGalerkin(GVec v, int n, GVec vecs[], int numFields, int *fields, NonlinearOperator op, PetscScalar alpha, PetscTruth isALE, void *ctx) { Grid grid; GVec x, y; int fieldIdx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidIntPointer(fields); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) { if ((fields[fieldIdx] < 0) || (fields[fieldIdx] > grid->numFields)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid field number"); } } switch(n) { case 1: x = y = vecs[0]; break; case 2: x = vecs[0]; y = vecs[1]; break; default: SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Only quadratic operators supported"); } /* Should check that fields are in the local order */ ierr = (*grid->ops->gvecevaluatenonlinearoperatorgalerkin)(grid, v, x, y, numFields, fields, grid->locOrder, op, alpha, isALE, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateOperatorGalerkin" /*@ GVecEvaluateOperatorGalerkin - Evaluates the weak form of an operator over a field on the locations defined by the underlying grid and its discretization. Input Parameter: . v - The grid vector . x,y - The input grid vectors, usually the previous solution . sField - The shape function field . tField - The test function field . op - The operator . alpha - A scalar multiplier . ctx - An optional user provided context for the function Level: intermediate .seealso: GVecEvaluateFunction,GMatEvaluateFunctionGalerkin @*/ int GVecEvaluateOperatorGalerkin(GVec v, GVec x, GVec y, int sField, int tField, int op, PetscScalar alpha, void *ctx) { Grid grid; VarOrdering sOldOrder, tOldOrder; VarOrdering sOrder, tOrder; LocalVarOrdering sLocOrder, tLocOrder; int numFields; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidHeaderSpecific(x, VEC_COOKIE); PetscValidHeaderSpecific(y, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); /* Check for valid fields */ ierr = GridGetNumFields(grid, &numFields); CHKERRQ(ierr); if ((sField < 0) || (sField >= numFields) || (tField < 0) || (tField >= numFields)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid field number"); } /* Check for valid operator */ if ((op < 0) || (op >= grid->fields[sField].disc->numOps) || (!grid->fields[sField].disc->operators[op])) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid operator"); } if ((grid->fields[sField].disc->operators[op]->precompInt == PETSC_NULL) && (grid->fields[sField].disc->operators[op]->opFunc == PETSC_NULL) && (grid->fields[sField].disc->operators[op]->ALEOpFunc == PETSC_NULL)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Operator cannot be calculated"); } if (grid->fields[sField].disc->numQuadPoints != grid->fields[tField].disc->numQuadPoints) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible quadrature sets"); } /* Create orderings */ ierr = GVecGetOrder(x, &sOldOrder); CHKERRQ(ierr); ierr = VarOrderingCreateSubset(sOldOrder, 1, &sField, PETSC_FALSE, &sOrder); CHKERRQ(ierr); ierr = LocalVarOrderingCreate(grid, 1, &sField, &sLocOrder); CHKERRQ(ierr); ierr = GVecGetOrder(v, &tOldOrder); CHKERRQ(ierr); ierr = VarOrderingCreateSubset(tOldOrder, 1, &tField, PETSC_FALSE, &tOrder); CHKERRQ(ierr); ierr = LocalVarOrderingCreate(grid, 1, &tField, &tLocOrder); CHKERRQ(ierr); /* Calculate operator */ ierr = (*grid->ops->gvecevaluateoperatorgalerkin)(grid, v, x, y, sOrder, sLocOrder, tOrder, tLocOrder, op, alpha, ctx); CHKERRQ(ierr); /* Destroy orderings */ ierr = VarOrderingDestroy(sOrder); CHKERRQ(ierr); ierr = LocalVarOrderingDestroy(sLocOrder); CHKERRQ(ierr); ierr = VarOrderingDestroy(tOrder); CHKERRQ(ierr); ierr = LocalVarOrderingDestroy(tLocOrder); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateOperatorGalerkinRectangular" /*@ GVecEvaluateOperatorGalerkinRectangular - Evaluates the weak form of an operator over a field on the locations defined by the underlying grid and its discretization. Collective on GVec Input Parameters: + v - The grid vector . x - A grid vector, usually the previous solution . sOrder - The global variable ordering for the shape functions . sLocOrder - The local variable ordering for the shape functions . sOrder - The global variable ordering for the test functions . sLocOrder - The local variable ordering for the test functions . op - The operator . alpha - The scalar multiplier for the operator - ctx - [Optional] The user provided context for the function Level: intermediate .keywords: grid vector, evaluate, operator, galerkin .seealso: GVecEvaluateFunction(), GMatEvaluateFunctionGalerkin() @*/ int GVecEvaluateOperatorGalerkinRectangular(GVec v, GVec x, VarOrdering sOrder, LocalVarOrdering sLocOrder, VarOrdering tOrder, LocalVarOrdering tLocOrder, int op, PetscScalar alpha, void *ctx) { Grid grid; FieldClassMap sMap, tMap; int numTotalFields; int *sFields, *tFields; int f; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidHeaderSpecific(x, VEC_COOKIE); PetscValidHeaderSpecific(sOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(sLocOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(tOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(tLocOrder, VAR_ORDER_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = GridGetNumFields(grid, &numTotalFields); CHKERRQ(ierr); ierr = VarOrderingGetClassMap(sOrder, &sMap); CHKERRQ(ierr); ierr = VarOrderingGetClassMap(tOrder, &tMap); CHKERRQ(ierr); sFields = sMap->fields; tFields = tMap->fields; /* Check for compatible orderings */ if ((tOrder->numVars != v->N) || (tOrder->numLocVars != v->n)) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function ordering"); } if ((sOrder->numVars != x->N) || (sOrder->numLocVars != x->n)) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible shape function ordering"); } if (sMap->numFields != tMap->numFields) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Shape and test function orderings must have the same number of fields"); } for(f = 0; f < sMap->numFields; f++) { if ((sFields[f] < 0) || (sFields[f] >= numTotalFields) || (tFields[f] < 0) || (tFields[f] >= numTotalFields)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid field number"); } /* Check for valid operator */ if ((op < 0) || (op >= grid->fields[sFields[f]].disc->numOps) || (!grid->fields[sFields[f]].disc->operators[op])) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid operator"); } if ((grid->fields[sFields[f]].disc->operators[op]->precompInt == PETSC_NULL) && (grid->fields[sFields[f]].disc->operators[op]->opFunc == PETSC_NULL) && (grid->fields[sFields[f]].disc->operators[op]->ALEOpFunc == PETSC_NULL)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid operator"); } if (grid->fields[sFields[f]].disc->numQuadPoints != grid->fields[tFields[f]].disc->numQuadPoints) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible quadrature sets"); } } /* Calculate operator */ ierr = (*grid->ops->gvecevaluateoperatorgalerkin)(grid, v, x, x, sOrder, sLocOrder, tOrder, tLocOrder, op, alpha, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateJacobian" /*@ GVecEvaluateJacobian - Evaluates the action of the system matrix over the active fields on the locations defined by the underlying grid and its discretization. Input Parameters: + x - The argument vector to the Jacobian . y - The vector to which the Jacobian is applied - ctx - An optional user provided context for the function Output Parameter: . v - The result vector J(x) y Level: intermediate .seealso: GVecEvaluateFunction, GMatEvaluateFunctionGalerkin @*/ int GVecEvaluateJacobian(GVec v, GVec x, GVec y, void *ctx) { Grid grid; Mesh mesh; MeshMover mover; Grid meshVelGrid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); 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); } ierr = PetscLogEventBegin(GVEC_EvaluateJacobian, v, grid, x, y); CHKERRQ(ierr); ierr = (*grid->ops->gvecevaluatesystemmatrix)(grid, x, y, v, ctx); CHKERRQ(ierr); ierr = PetscLogEventEnd(GVEC_EvaluateJacobian, v, grid, x, y); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateJacobianDiagonal" /*@ GVecEvaluateJacobianDiagonal - Evaluates the action of the diagonal of the system matrix over the active fields on the locations defined by the underlying grid and its discretization. Input Parameters: + x - The argument vector to the Jacobian - ctx - An optional user provided context for the function Output Parameter: . d - The diagonal of J(x) Level: intermediate .seealso: GVecEvaluateJacobian, GVecEvaluateFunction @*/ int GVecEvaluateJacobianDiagonal(GVec d, GVec x, void *ctx) { Grid grid; Mesh mesh; MeshMover mover; Grid meshVelGrid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(d, VEC_COOKIE); ierr = GVecGetGrid(d, &grid); CHKERRQ(ierr); 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); } ierr = PetscLogEventBegin(GVEC_EvaluateJacobian, d, grid, x, 0); CHKERRQ(ierr); ierr = (*grid->ops->gvecevaluatesystemmatrixdiagonal)(grid, x, d, ctx); CHKERRQ(ierr); ierr = PetscLogEventEnd(GVEC_EvaluateJacobian, d, grid, x, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecEvaluateJacobianConstrained" /*@ GVecEvaluateJacobianConstrained - Evaluates the action of the Jacobian of the extra fields introduced by the constraints. Collective on GVec Input Parameter: . v - The grid vector Output Parameter . x - The action of the constraint Jacobian Level: intermediate .keywords: grid vector, jacobian, evaluate, constraint .seealso: GVecEvaluateFunction(), GVecEvaluateRhs() @*/ int GVecEvaluateJacobianConstrained(GVec v, GVec x) { Grid grid; PetscTruth isConstrained; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidHeaderSpecific(x, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = GridIsConstrained(grid, &isConstrained); CHKERRQ(ierr); if ((isConstrained == PETSC_TRUE) && (grid->constraintCtx->ops->applyjac)) { ierr = (*grid->constraintCtx->ops->applyjac)(grid, x, v); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecSolveJacobianConstrained" /*@ GVecSolveJacobianConstrained - Evaluates the action of the inverse of the Jacobian of the extra fields introduced by the constraints. Collective on GVec Input Parameter: . v - The grid vector Output Parameter: . x - The action of the inverse of the constraint Jacobian Level: intermediate .keywords: grid vector, jacobian, solve, constraint .seealso: GVecEvaluateFunction(), GVecEvaluateRhs() @*/ int GVecSolveJacobianConstrained(GVec v, GVec x) { Grid grid; PetscTruth isConstrained; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidHeaderSpecific(x, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = GridIsConstrained(grid, &isConstrained); CHKERRQ(ierr); if ((isConstrained == PETSC_TRUE) && (grid->constraintCtx->ops->solvejac)) { ierr = (*grid->constraintCtx->ops->solvejac)(grid, x, v); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecSetBoundary" /*@ GVecSetBoundary - Applies the specified Dirichlet boundary conditions to this vector. Input Parameter: . v - The grid vector . ctx - An optional user provided context for the function Level: intermediate .seealso: GridSetBC(), GridAddBC() @*/ int GVecSetBoundary(GVec v, void *ctx) { Grid grid; int bc; int num; int *bd; int *field; PointFunction *func; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = PetscLogEventBegin(GVEC_SetBoundary, v, grid, 0, 0); CHKERRQ(ierr); for(bc = 0, num = 0; bc < grid->numBC; bc++) { if (grid->bc[bc].reduce == PETSC_FALSE) num++; } if (num > 0) { ierr = PetscMalloc(num * sizeof(int), &bd); CHKERRQ(ierr); ierr = PetscMalloc(num * sizeof(int), &field); CHKERRQ(ierr); ierr = PetscMalloc(num * sizeof(PointFunction *), &func); CHKERRQ(ierr); for(bc = 0, num = 0; bc < grid->numBC; bc++) { if (grid->bc[bc].reduce == PETSC_FALSE) { bd[num] = grid->bc[bc].boundary; field[num] = grid->bc[bc].field; func[num++] = grid->bc[bc].func; } } ierr = GridSetVecBoundaryRectangular(num, bd, field, func, grid->order, v, ctx); CHKERRQ(ierr); ierr = PetscFree(bd); CHKERRQ(ierr); ierr = PetscFree(field); CHKERRQ(ierr); ierr = PetscFree(func); CHKERRQ(ierr); } for(bc = 0; bc < grid->numPointBC; bc++) { if (grid->pointBC[bc].reduce == PETSC_FALSE) { ierr = GridSetVecPointBoundary(grid->pointBC[bc].node, grid->pointBC[bc].field, grid->pointBC[bc].func, v, ctx);CHKERRQ(ierr); } } ierr = PetscLogEventEnd(GVEC_SetBoundary, v, grid, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecSetBoundaryZero" /*@ GVecSetBoundaryZero - Applies a Dirichlet boundary condition of zero to the specified fields in this vector. This is mainly used to double up the boundary conditions in a time dependent nonlinear problem. The original boundary condition is used to initialize the solution at the start of the Newton iteration, and this is used to set the corresponding components of the residual to zero. Input Parameter: . v - The grid vector . ctx - An optional user provided context for the function Level: intermediate .seealso: GridSetBC(), GridAddBC() @*/ int GVecSetBoundaryZero(GVec v, void *ctx) { Grid grid; int bc; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); for(bc = 0; bc < grid->numBC; bc++) { if (grid->bc[bc].reduce == PETSC_FALSE) { ierr = GridSetVecBoundary(grid->bc[bc].boundary, grid->bc[bc].field, PointFunctionZero, v, ctx); CHKERRQ(ierr); } } for(bc = 0; bc < grid->numPointBC; bc++) { if (grid->pointBC[bc].reduce == PETSC_FALSE) { ierr = GridSetVecPointBoundary(grid->pointBC[bc].node, grid->pointBC[bc].field, PointFunctionZero, v, ctx);CHKERRQ(ierr); } } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecSetBoundaryDifference" /*@ GVecSetBoundaryDifference - Applies the specified Dirichlet boundary conditions to this vector, but uses the difference of the value in the given vector and the computed value. Input Parameter: . v - The grid vector . u - A grid vector, usually the previous solution . ctx - An optional user provided context for the function Level: intermediate .seealso: GridSetVecBoundaryDifference(), GridSetBC(), GridAddBC() @*/ int GVecSetBoundaryDifference(GVec v, GVec u, void *ctx) { Grid grid; int bc; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); PetscValidHeaderSpecific(u, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); for(bc = 0; bc < grid->numBC; bc++) { if (grid->bc[bc].reduce == PETSC_FALSE) { ierr = GridSetVecBoundaryDifference(grid->bc[bc].boundary, grid->bc[bc].field, u, grid->bc[bc].func, v, ctx);CHKERRQ(ierr); } } for(bc = 0; bc < grid->numPointBC; bc++) { if (grid->pointBC[bc].reduce == PETSC_FALSE) { ierr = GridSetVecPointBoundaryDifference(grid->pointBC[bc].node, grid->pointBC[bc].field, u, grid->pointBC[bc].func, v, ctx);CHKERRQ(ierr); } } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecInterpolateField" /*@ GVecInterpolateField - Interpolates the field defined by the grid vector to the specified point, and stores the components in the array provided. Input Parameter: . v - The grid vector . x,y,z - The interpolation point . ctx - An InterpolationContext Output Parameter: . values - The interpolated field at the point Level: intermediate .keywords: grid vector, interpolation .seealso: GVecEvaluateFunction, GVecEvaluateFunctionGalerkin @*/ int GVecInterpolateField(GVec v, int field, double x, double y, double z, PetscScalar *values, InterpolationContext *ctx) { Grid grid; Discretization disc; int *elemStart; ElementVec vec; PetscScalar *array; double coords[3]; double *tmp; MPI_Status status; PetscTruth lost; /*double lostDist, lostX, lostY, lostZ;*/ int numFound, numSent; int comp, elem, newElem, closeNode; int proc, var; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(v, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = PetscLogEventBegin(GVEC_InterpolateField, v, grid, 0, 0); CHKERRQ(ierr); disc = grid->fields[field].disc; comp = disc->comp; elemStart = grid->locOrder->elemStart; vec = grid->ghostElementVec; array = &vec->array[elemStart[field]]; /* Watch out for null collective call */ if (values == PETSC_NULL) { /* Mark point as inactive */ elem = -2; } else { /* Locate the point in the mesh */ ierr = MeshLocatePoint(ctx->mesh, x, y, 0.0, &elem); CHKERRQ(ierr); } if (elem >= 0) { /* Initialize element vector */ ierr = ElementVecZero(vec); CHKERRQ(ierr); vec->reduceSize = grid->locOrder->elemSize; /* Get indices for the old element */ ierr = GridCalcInterpolationElementVecIndices(grid, ctx->mesh, elem, ctx->order, vec); CHKERRQ(ierr); /* Get values from local vector to element vector */ ierr = GridLocalToElementGeneral(grid, ctx->ghostVec, grid->bdReduceVecCur, grid->reduceSystem, grid->reduceElement, vec);CHKERRQ(ierr); if (grid->explicitConstraints == PETSC_TRUE) { /* Must transform to unconstrained variables for element integrals */ ierr = GridProjectInterpolationElementVec(grid, ctx->mesh, elem, ctx->order, PETSC_FALSE, vec); CHKERRQ(ierr); } /* Interpolate data to get field at new point */ ierr = DiscretizationInterpolateField(disc, ctx->mesh, elem, x, y, 0.0, array, values, INTERPOLATION_LOCAL); CHKERRQ(ierr); } else if (elem != -2) { if (ctx->numProcs == 1) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Node (%g,%g) not found in new mesh", x, y); } else if (ctx->batchMode == PETSC_TRUE) { if (ctx->numPoints[ctx->rank] >= ctx->maxPoints) { ierr = PetscMalloc(ctx->maxPoints*6 * sizeof(double), &tmp); CHKERRQ(ierr); ierr = PetscMemcpy(tmp, ctx->coords, ctx->maxPoints*3 * sizeof(double)); CHKERRQ(ierr); ierr = PetscFree(ctx->coords); CHKERRQ(ierr); ctx->maxPoints *= 2; ctx->coords = tmp; } ctx->coords[ctx->numPoints[ctx->rank]*3] = x; ctx->coords[ctx->numPoints[ctx->rank]*3+1] = y; ctx->coords[ctx->numPoints[ctx->rank]*3+2] = z; ctx->numPoints[ctx->rank]++; for(var = 0; var < comp; var++) values[var] = 0.0; } } if ((ctx->numProcs > 1) && (ctx->batchMode == PETSC_FALSE)) { /* Communicate missing nodes */ ierr = MPI_Allgather(&elem, 1, MPI_INT, ctx->activeProcs, 1, MPI_INT, grid->comm); CHKERRQ(ierr); coords[0] = x; coords[1] = y; coords[2] = z; ierr = MPI_Allgather(coords, 3, MPI_DOUBLE, ctx->coords, 3, MPI_DOUBLE, grid->comm); CHKERRQ(ierr); /* Handle each point */ numFound = 0; for(proc = 0; proc < ctx->numProcs; proc++) { /* If the point was not located by someone else */ if ((ctx->activeProcs[proc] == -1) && (proc != ctx->rank)) { ierr = MeshLocatePoint(ctx->mesh, ctx->coords[proc*3], ctx->coords[proc*3+1], 0.0, &newElem); CHKERRQ(ierr); if (newElem >= 0) { /* Mark point as found by this processor */ ctx->activeProcs[proc] = ctx->rank; /* Get indices for the old element */ ierr = GridCalcInterpolationElementVecIndices(grid, ctx->mesh, newElem, ctx->order, vec); CHKERRQ(ierr); /* Get values from local vector to element vector */ ierr = GridLocalToElementGeneral(grid, ctx->ghostVec, grid->bdReduceVecCur, grid->reduceSystem, grid->reduceElement, vec);CHKERRQ(ierr); if (grid->explicitConstraints == PETSC_TRUE) { /* Must transform to unconstrained variables for element integrals */ ierr = GridProjectInterpolationElementVec(grid, ctx->mesh, newElem, ctx->order, PETSC_FALSE, vec); CHKERRQ(ierr); } /* Interpolate data to get field at new point */ ierr = DiscretizationInterpolateField(disc, ctx->mesh, newElem, ctx->coords[proc*3], ctx->coords[proc*3+1], 0.0, array, &ctx->values[proc*comp], INTERPOLATION_LOCAL); CHKERRQ(ierr); numFound++; } } else { /* Mark point as inactive */ ctx->activeProcs[proc] = -2; } } /* Extra step to handle nodes which slip outside the boundary */ ierr = MPI_Allreduce(ctx->activeProcs, ctx->calcProcs, ctx->numProcs, MPI_INT, MPI_MAX, grid->comm); CHKERRQ(ierr); lost = PETSC_FALSE; for(proc = 0; proc < ctx->numProcs; proc++) { if (ctx->calcProcs[proc] == -1) { /* PetscPrintf(grid->comm, "Lost point:\n "); for(var = 0; var < ctx->numProcs; var++) PetscPrintf(grid->comm, " aProc[%d]: %d", var, ctx->activeProcs[var]); PetscPrintf(grid->comm, "\n"); for(var = 0; var < ctx->numProcs; var++) PetscPrintf(grid->comm, " cProc[%d]: %d", var, ctx->calcProcs[var]); PetscPrintf(grid->comm, "\n"); PetscPrintf(grid->comm, " (%g, %g)\n", ctx->coords[proc*3], ctx->coords[proc*3+1]);*/ lost = PETSC_TRUE; /* No one found it so search for the closest element */ ierr = MeshGetNearestNode(ctx->mesh, ctx->coords[proc*3], ctx->coords[proc*3+1], 0.0, PETSC_TRUE, &closeNode); CHKERRQ(ierr); if (closeNode >= 0) { ierr = MeshGetElementFromNode(ctx->mesh, closeNode, &newElem); CHKERRQ(ierr); /* PetscPrintf(PETSC_COMM_SELF, " found lost point on processor %d in element %d\n", ctx->rank, newElem); ierr = MeshGetNodeCoords(ctx->mesh, closeNode, &lostX, &lostY, &lostZ); CHKERRQ(ierr); lostDist = PetscSqr(MeshPeriodicDiffX(ctx->mesh, lostX - ctx->coords[proc*3])) + PetscSqr(MeshPeriodicDiffY(ctx->mesh, lostY - ctx->coords[proc*3+1])); PetscPrintf(PETSC_COMM_SELF, " at (%g, %g), it was %g away from the point\n", lostX, lostY, lostDist); */ /* Mark point as found by this processor */ ctx->activeProcs[proc] = ctx->rank; /* Get indices for the old element */ ierr = GridCalcInterpolationElementVecIndices(grid, ctx->mesh, newElem, ctx->order, vec); CHKERRQ(ierr); /* Get values from local vector to element vector */ ierr = GridLocalToElementGeneral(grid, ctx->ghostVec, grid->bdReduceVecCur, grid->reduceSystem, grid->reduceElement, vec); CHKERRQ(ierr); if (grid->explicitConstraints == PETSC_TRUE) { /* Must transform to unconstrained variables for element integrals */ ierr = GridProjectInterpolationElementVec(grid, ctx->mesh, newElem, ctx->order, PETSC_FALSE, vec); CHKERRQ(ierr); } /* Interpolate data to get field at new point */ ierr = DiscretizationInterpolateField(disc, ctx->mesh, newElem, ctx->coords[proc*3], ctx->coords[proc*3+1], 0.0, array, &ctx->values[proc*comp], INTERPOLATION_LOCAL); CHKERRQ(ierr); numFound++; } } } /* Communicate interpolated values */ if (lost == PETSC_TRUE) { ierr = MPI_Allreduce(ctx->activeProcs, ctx->calcProcs, ctx->numProcs, MPI_INT, MPI_MAX, grid->comm); CHKERRQ(ierr); /* PetscPrintf(grid->comm, "Recommunicated interpolated values\n "); for(var = 0; var < ctx->numProcs; var++) PetscPrintf(grid->comm, " aProc[%d]: %d", var, ctx->activeProcs[var]); PetscPrintf(grid->comm, "\n"); for(var = 0; var < ctx->numProcs; var++) PetscPrintf(grid->comm, " cProc[%d]: %d", var, ctx->calcProcs[var]); PetscPrintf(grid->comm, "\n"); */ } numSent = 0; for(proc = 0; proc < ctx->numProcs; proc++) { if ((elem == -1) && (ctx->rank == proc)) { /* Point from this domain was not found */ if (ctx->calcProcs[proc] < 0) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Node not found at (%g,%g)", x, y); } /* Must check for lost points found in the same domain */ if (ctx->calcProcs[proc] != ctx->rank) { ierr = MPI_Recv(values, comp, MPI_DOUBLE, ctx->calcProcs[proc], 0, grid->comm, &status); CHKERRQ(ierr); } else { for(var = 0; var < comp; var++) values[var] = ctx->values[proc*comp+var]; numSent++; } #if 0 PetscPrintf(PETSC_COMM_SELF, "Immediate point %d on proc %d\n", ctx->curPoint++, ctx->rank); PetscPrintf(PETSC_COMM_SELF, " x: %g y: %g z: %g\n ", ctx->coords[proc*3], ctx->coords[proc*3+1], ctx->coords[proc*3+2]); for(c = 0; c < comp; c++) { PetscPrintf(PETSC_COMM_SELF, "val[%d]: %g ", c, values[c]); } PetscPrintf(PETSC_COMM_SELF, "\n"); #endif } else if (ctx->rank == ctx->calcProcs[proc]) { /* Point from another domain was found here */ ierr = MPI_Send(&ctx->values[proc*comp], comp, MPI_DOUBLE, proc, 0, grid->comm); CHKERRQ(ierr); numSent++; } else if (ctx->rank == ctx->activeProcs[proc]) { /* Point was found by multiple processors */ if (ctx->calcProcs[proc] < ctx->rank) { SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid precendence during interpolation"); } numSent++; } } if (numSent != numFound) { SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Number of nodes found %d does not match number sent %d", numFound, numSent); } /* if (lost == PETSC_TRUE) { PetscBarrier((PetscObject) grid); PetscPrintf(grid->comm, "Finished interpolation after lost point\n"); }*/ } ierr = PetscLogEventEnd(GVEC_InterpolateField, v, grid, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecInterpolateFieldBatch" int GVecInterpolateFieldBatch(GVec v, int field, InterpolationContext *ctx) { Grid grid; Discretization disc; int comp; ElementVec vec; int numProcs = ctx->numProcs; int rank = ctx->rank; int *numCoords, *cOffsets; int *activePoints, *calcPoints; double *coords; PetscScalar *values, *array; int totPoints; int numFound; int proc, p, elem; int ierr; PetscFunctionBegin; if ((ctx->numProcs == 1) || (ctx->batchMode == PETSC_FALSE)) PetscFunctionReturn(0); PetscValidHeaderSpecific(v, VEC_COOKIE); ierr = GVecGetGrid(v, &grid); CHKERRQ(ierr); ierr = PetscLogEventBegin(GVEC_InterpolateFieldBatch, v, grid, 0, 0); CHKERRQ(ierr); disc = grid->fields[field].disc; comp = disc->comp; vec = grid->ghostElementVec; array = &vec->array[grid->locOrder->elemStart[field]]; /* Get the number of points from each domain */ ierr = PetscMalloc((numProcs+1) * sizeof(int), &numCoords); CHKERRQ(ierr); ierr = PetscMalloc((numProcs+1) * sizeof(int), &cOffsets); CHKERRQ(ierr); ierr = PetscLogEventBegin(GVEC_InterpolateFieldBatchParallel, v, grid, 0, 0); CHKERRQ(ierr); ierr = MPI_Allgather(&ctx->numPoints[rank], 1, MPI_INT, &ctx->numPoints[1], 1, MPI_INT, v->comm); CHKERRQ(ierr); ierr = PetscLogEventEnd(GVEC_InterpolateFieldBatchParallel, v, grid, 0, 0); CHKERRQ(ierr); /* Get the coordintes of each point */ ctx->numPoints[0] = 0; cOffsets[0] = 0; for(proc = 0; proc < numProcs; proc++) { numCoords[proc] = ctx->numPoints[proc+1]*3; ctx->numPoints[proc+1] += ctx->numPoints[proc]; cOffsets[proc+1] = ctx->numPoints[proc+1]*3; } totPoints = ctx->numPoints[numProcs]; ierr = PetscMalloc(totPoints*3 * sizeof(double), &coords); CHKERRQ(ierr); ierr = PetscLogEventBegin(GVEC_InterpolateFieldBatchParallel, v, grid, 0, 0); CHKERRQ(ierr); ierr = MPI_Allgatherv(ctx->coords, numCoords[rank], MPI_DOUBLE, coords, numCoords, cOffsets, MPI_DOUBLE, v->comm); CHKERRQ(ierr); ierr = PetscLogEventEnd(GVEC_InterpolateFieldBatchParallel, v, grid, 0, 0); CHKERRQ(ierr); ierr = PetscMalloc(totPoints * sizeof(int), &activePoints); CHKERRQ(ierr); ierr = PetscMalloc(totPoints * sizeof(int), &calcPoints); CHKERRQ(ierr); ierr = PetscMalloc(totPoints*comp * sizeof(PetscScalar), &values); CHKERRQ(ierr); ierr = PetscMalloc(totPoints*comp * sizeof(PetscScalar), &ctx->values); CHKERRQ(ierr); ierr = PetscMemzero(values, totPoints*comp * sizeof(PetscScalar)); CHKERRQ(ierr); ierr = PetscLogEventBegin(GVEC_InterpolateFieldBatchCalc, v, grid, 0, 0); CHKERRQ(ierr); /* Mark this domain's points as inactive */ for(p = 0; p < totPoints; p++) { if ((p >= ctx->numPoints[rank]) && (p < ctx->numPoints[rank+1])) activePoints[p] = -2; else activePoints[p] = -1; } /* Handle each point */ for(p = 0, numFound = 0; p < totPoints; p++) { if (activePoints[p] > -2) { ierr = MeshLocatePoint(ctx->mesh, coords[p*3], coords[p*3+1], coords[p*3+2], &elem); CHKERRQ(ierr); if (elem >= 0) { /* Mark point as found by this processor */ activePoints[p] = rank; /* Get indices for the old element */ ierr = GridCalcInterpolationElementVecIndices(grid, ctx->mesh, elem, ctx->order, vec); CHKERRQ(ierr); /* Get values from local vector to element vector */ ierr = GridLocalToElementGeneral(grid, ctx->ghostVec, grid->bdReduceVecCur, grid->reduceSystem, grid->reduceElement, vec); CHKERRQ(ierr); if (grid->explicitConstraints == PETSC_TRUE) { /* Must transform to unconstrained variables for element integrals */ ierr = GridProjectInterpolationElementVec(grid, ctx->mesh, elem, ctx->order, PETSC_FALSE, vec); CHKERRQ(ierr); } /* Interpolate data to get field at new point */ ierr = DiscretizationInterpolateField(disc, ctx->mesh, elem, coords[p*3], coords[p*3+1], coords[p*3+2], array, &values[p*comp], INTERPOLATION_LOCAL); CHKERRQ(ierr); numFound++; } } } ierr = PetscLogEventEnd(GVEC_InterpolateFieldBatchCalc, v, grid, 0, 0); CHKERRQ(ierr); /* Communicate interpolated values */ #if 1 /* This is REALLY bad but probably better than what is there now */ ierr = PetscLogEventBegin(GVEC_InterpolateFieldBatchParallel, v, grid, 0, 0); CHKERRQ(ierr); ierr = MPI_Allreduce(values, ctx->values, totPoints*comp, MPI_DOUBLE, MPI_SUM, v->comm); CHKERRQ(ierr); ierr = PetscLogEventEnd(GVEC_InterpolateFieldBatchParallel, v, grid, 0, 0); CHKERRQ(ierr); #else ierr = MPI_Allreduce(activePoints, calcPoints, totPoints, MPI_INT, MPI_MAX, v->comm); CHKERRQ(ierr); for(proc = 0, numSent = 0; proc < numProcs; proc++) for(p = ctx->numPoints[proc]; p < ctx->numPoints[proc+1]; p++) { if (rank == proc) { /* Point from this domain was not found */ if (calcPoints[p] < 0) { PetscLogInfo(v, "[%d]x: %g y: %g z: %g\n", rank, coords[p*3], coords[p*3+1], coords[p*3+2]); SETERRQ(PETSC_ERR_PLIB, "Node not found"); } ierr = MPI_Recv(&ctx->values[p*comp], comp, MPI_DOUBLE, calcPoints[p], 0, v->comm, &status); CHKERRQ(ierr); } else if (rank == calcPoints[p]) { /* Point from another domain was found here */ ierr = MPI_Send(&values[p*comp], comp, MPI_DOUBLE, proc, 0, v->comm); CHKERRQ(ierr); numSent++; } else if (rank == activePoints[p]) { /* Point was found by multiple processors */ if (calcPoints[p] < rank) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid precendence during interpolation"); numSent++; } } if (numSent != numFound) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Number of nodes found does not match number sent"); #endif ierr = PetscFree(numCoords); CHKERRQ(ierr); ierr = PetscFree(activePoints); CHKERRQ(ierr); ierr = PetscFree(calcPoints); CHKERRQ(ierr); ierr = PetscFree(values); CHKERRQ(ierr); ierr = PetscLogEventEnd(GVEC_InterpolateFieldBatch, v, grid, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecCreate" /*@ GVecCreate - Creates a grid vector associated with a particular grid. Input Parameter: . grid - The grid defining the discrete function Output Parameter: . gvec - The resulting grid vector Level: beginner .seealso GVecCreateConstrained @*/ int GVecCreate(Grid grid, GVec *gvec) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(gvec); ierr = GridSetUp(grid); if (ierr) SETERRQ(PETSC_ERR_PLIB, "Grid setup failed"); ierr = GVecCreateRectangular(grid, grid->order, gvec); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecCreateGhost" /*@ GVecCreateGhost - Creates a grid vector associated with a particular grid with ghost padding on each processor. Input Parameter: . grid - The grid defining the discrete function Output Parameter: . gvec - The resulting grid vector Level: beginner .seealso GVecCreate, GVecCreateConstrained @*/ int GVecCreateGhost(Grid grid, GVec *gvec) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(gvec); ierr = GridSetUp(grid); if (ierr) SETERRQ(PETSC_ERR_PLIB, "Grid setup failed"); ierr = GVecCreateRectangularGhost(grid, grid->order, gvec); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecCreateRectangular" /*@ GVecCreateRectangular - Creates a grid vector associated with a particular grid and ordering. Input Parameter: . grid - The grid defining the discrete function . order - The variable ordering Output Parameter: . gvec - The resulting grid vector Level: beginner .seealso GVecCreate, GVecCreateConstrained @*/ int GVecCreateRectangular(Grid grid, VarOrdering order, GVec *gvec) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE); PetscValidPointer(gvec); ierr = GridSetUp(grid); if (ierr) SETERRQ(PETSC_ERR_PLIB, "Grid setup failed"); ierr = VecCreateMPI(grid->comm, order->numLocVars, order->numVars, gvec); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *gvec, "Grid", (PetscObject) grid); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *gvec, "Order", (PetscObject) order); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecCreateRectangularGhost" /*@ GVecCreateRectangularGhost - Creates a grid vector associated with a particular grid and ordering, adding ghost padding on each processor. Input Parameter: . grid - The grid defining the discrete function . order - The variable ordering Output Parameter: . gvec - The resulting grid vector Level: beginner .seealso GVecCreate, GVecCreateConstrained @*/ int GVecCreateRectangularGhost(Grid grid, VarOrdering order, GVec *gvec) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE); PetscValidPointer(gvec); ierr = GridSetUp(grid); if (ierr) SETERRQ(PETSC_ERR_PLIB, "Grid setup failed"); /* I should use MPICreateGhost(), but it is in a state of flux */ ierr = VecCreate(grid->comm, gvec); CHKERRQ(ierr); ierr = VecSetSizes(*gvec, order->numLocVars, order->numVars); CHKERRQ(ierr); ierr = VecCreate_MPI_Private(*gvec, order->numOverlapVars-order->numLocVars, PETSC_NULL, PETSC_NULL); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *gvec, "Grid", (PetscObject) grid); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *gvec, "Order", (PetscObject) order); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecCreateConstrained" /*@ GVecCreateConstrained - Creates a grid vector associated with a particular grid, but after constraints have been applied Input Parameter: . grid - The grid defining the discrete function Output Parameter: . gvec - The resulting grid vector Level: beginner .seealso GVecCreate @*/ int GVecCreateConstrained(Grid grid, GVec *gvec) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(gvec); ierr = GridSetUp(grid); if (ierr) SETERRQ(PETSC_ERR_PLIB, "Grid setup failed"); if (grid->isConstrained) { ierr = GVecCreateRectangular(grid, grid->constraintOrder, gvec); CHKERRQ(ierr); } else { ierr = GVecCreate(grid, gvec); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecCreateBoundaryRestriction" /*@ GVecCreateBoundaryRestriction - Creates a grid vector associated with the boundary variables of a particular grid. Input Parameter: . grid - The grid defining the discrete function Output Parameter: . gvec - The resulting grid vector Level: beginner .seealso GVecCreate, GVecCreateConstrained @*/ int GVecCreateBoundaryRestriction(Grid grid, GVec *gvec) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(gvec); ierr = GridSetUp(grid); if (ierr) SETERRQ(PETSC_ERR_PLIB, "Grid setup failed"); ierr = GridSetupBoundary(grid); CHKERRQ(ierr); ierr = GVecCreateRectangular(grid, grid->bdOrder, gvec); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "PointFunctionInterpolateField" /*@ PointFunctionInterpolateField - A PointFunction that gives the field value at a given point in the mesh based upon a grid vector representation. Input Paramters: . n - The number of points . comp - The number of components . x,y,z - The coordinates of points . ctx - An InterpolationContext Output Paramter: . values - The location where the field values are stored Level: advanced .keywords point function, grid .seealso PointFunctionOne, PointFunctionZero, PointFunctionConstant @*/ int PointFunctionInterpolateField(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) { InterpolationContext *iCtx = (InterpolationContext *) ctx; double *px = x; double *py = y; double *pz = z; int i; int ierr; PetscFunctionBegin; /* Check for a call purely for collective operations */ if (n == 0) { ierr = GVecInterpolateField(iCtx->vec, iCtx->field, 0.0, 0.0, 0.0, PETSC_NULL, iCtx); CHKERRQ(ierr); PetscFunctionReturn(0); } /* We might do better by assuming that any two of these could be null */ if (z == PETSC_NULL) pz = px; if (y == PETSC_NULL) py = px; for(i = 0; i < n; i++) { ierr = GVecInterpolateField(iCtx->vec, iCtx->field, px[i], py[i], pz[i], &values[i*comp], iCtx); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "PointFunctionInterpolateFieldBatch" /*@ PointFunctionInterpolateFieldBatch - A PointFunction that finishes an interpolation by setting the off-processor values stored in the InterpolationContext. Input Paramters: . n - The number of points . comp - The number of components . x,y,z - The coordinates of points . ctx - An InterpolationContext Output Paramter: . values - The location where the field values are stored Level: advanced .keywords point function, grid .seealso PointFunctionOne, PointFunctionZero, PointFunctionConstant @*/ int PointFunctionInterpolateFieldBatch(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) { InterpolationContext *iCtx = (InterpolationContext *) ctx; double *px = x; double *py = y; double *pz = z; int p = iCtx->curPoint - iCtx->numPoints[iCtx->rank]; int i, c; PetscFunctionBegin; /* We might do better by assuming that any two of these could be null */ if (z == PETSC_NULL) pz = px; if (y == PETSC_NULL) py = px; for(i = 0; i < n; i++) { if ((*px == iCtx->coords[p*3]) && (*py == iCtx->coords[p*3+1]) && (*pz == iCtx->coords[p*3+2])) { /* Copy stored values */ #if 0 PetscPrintf(PETSC_COMM_SELF, "Batch point %d on proc %d\n", iCtx->curPoint, iCtx->rank); PetscPrintf(PETSC_COMM_SELF, " x: %g y: %g z: %g\n ", px[i], py[i], pz[i]); for(c = 0; c < comp; c++) { values[i*comp+c] = iCtx->values[iCtx->curPoint*comp+c]; PetscPrintf(PETSC_COMM_SELF, "val[%d]: %g ", c, values[i*comp+c]); } PetscPrintf(PETSC_COMM_SELF, "\n"); #else for(c = 0; c < comp; c++) values[i*comp+c] = iCtx->values[iCtx->curPoint*comp+c]; #endif iCtx->curPoint++; p++; } } PetscFunctionReturn(0); }