#ifdef PETSC_RCS_HEADER static char vcid[] = "$Id: gridBC.c,v 1.3 2000/01/30 18:27:13 huangp Exp $"; #endif #include "src/grid/gridimpl.h" /*I "grid.h" I*//*I "gvec.h" I*/ #include "src/vec/impls/mpi/pvecimpl.h" /* For GridCalcBCValues(), soon will not be needed */ /*------------------------------------------------ Standard Functions -----------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridDuplicateBC" /*@ GridDuplicateBC - Duplicates the boundary conditions of one grid in another. Collective on Grid Input Parameter: . grid - The grid Output Parameter: . newGrid - The altered grid Level: intermediate .keywords: grid, duplicate, BC .seealso: GridDuplicate() @*/ int GridDuplicateBC(Grid grid, Grid newGrid) { int bc; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidHeaderSpecific(newGrid, GRID_COOKIE); for(bc = 0; bc < grid->numBC; bc++) { ierr = GridAddBC(newGrid, grid->bc[bc].boundary, grid->bc[bc].field, grid->bc[bc].func, grid->bc[bc].reduce);CHKERRQ(ierr); } for(bc = 0; bc < grid->numPointBC; bc++) { ierr = GridAddPointBC(newGrid, grid->bc[bc].point[0], grid->bc[bc].point[1], grid->bc[bc].point[2], grid->bc[bc].field, grid->bc[bc].func, grid->bc[bc].reduce);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridFinalizeBC" /*@ GridFinalizeBC - Destroys all structures associated with explicit system reduction using boundary conditions. This should be called after all calculation is finished, prior to GridDestroy(). Collective on Grid Input Parameter: . grid - The grid Level: beginner .keywords: grid, boundary conditions .seealso: GridSetBC(), GridAddBC() @*/ int GridFinalizeBC(Grid grid) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); if (grid->bdReduceVec) { ierr = GVecDestroy(grid->bdReduceVec); CHKERRQ(ierr); } if (grid->bdReduceVecOld) { ierr = GVecDestroy(grid->bdReduceVecOld); CHKERRQ(ierr); } if (grid->bdReduceVecDiff) { ierr = GVecDestroy(grid->bdReduceVecDiff); CHKERRQ(ierr); } if (grid->bdReduceMat) { ierr = GMatDestroy(grid->bdReduceMat); CHKERRQ(ierr); } if (grid->reduceVec) { ierr = GVecDestroy(grid->reduceVec); CHKERRQ(ierr); } PetscFunctionReturn(0); } /*----------------------------------------------- Database Functions ------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridSetBC" /*@C GridSetBC This function sets the boundary condition to use for the problem. Collective on Grid Input Parameters: + grid - The grid . bd - The marker for the boundary along which conditions are applied . field - The field to which the boundary condition is applied . f - The function which defines the boundary condition - reduce - The flag for explicit reduction of the system Level: intermediate .keywords active field .seealso GridAddBC, GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction @*/ int GridSetBC(Grid grid, int bd, int field, PointFunction f, PetscTruth reduce) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); GridValidField(grid, field); grid->numBC = 0; ierr = GridAddBC(grid, bd, field, f, reduce); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridAddBC" /*@C GridAddBC This function adds a boundary condition to use for the problem. Collective on Grid Input Parameters: + grid - The grid . bd - The marker for the boundary along which conditions are applied . field - The field to which the boundary condition is applied . f - The function which defines the boundary condition - reduce - The flag for explicit reduction of the system Level: intermediate .keywords active field .seealso GridSetBC, GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction @*/ int GridAddBC(Grid grid, int bd, int field, PointFunction f, PetscTruth reduce) { GridBC *tempBC; int bdIndex; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); GridValidField(grid, field); while (grid->numBC >= grid->maxBC) { ierr = PetscMalloc(grid->maxBC*2 * sizeof(GridBC), &tempBC); CHKERRQ(ierr); PetscLogObjectMemory(grid, grid->maxBC * sizeof(GridBC)); ierr = PetscMemcpy(tempBC, grid->bc, grid->maxBC * sizeof(GridBC)); CHKERRQ(ierr); ierr = PetscFree(grid->bc); CHKERRQ(ierr); grid->bc = tempBC; grid->maxBC *= 2; } /* Make sure boundary is legal */ ierr = MeshGetBoundaryIndex(grid->mesh, bd, &bdIndex); CHKERRQ(ierr); grid->bc[grid->numBC].boundary = bd; grid->bc[grid->numBC].field = field; grid->bc[grid->numBC].func = f; grid->bc[grid->numBC].reduce = reduce; grid->bc[grid->numBC].node = -1; grid->numBC++; /* Check whether to reduce system */ if (reduce == PETSC_TRUE) grid->reduceSystem = PETSC_TRUE; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetPointBC" /*@C GridSetPointBC This function sets the boundary condition to use for the problem at a point. Collective on Grid Input Parameters: + grid - The grid . x,y,z - The point at which conditions are applied . field - The field to which the boundary condition is applied . f - The function which defines the boundary condition - reduce - The flag for explicit reduction of the system Level: intermediate .keywords active field .seealso GridAddBC, GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction @*/ int GridSetPointBC(Grid grid, double x, double y, double z, int field, PointFunction f, PetscTruth reduce) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); grid->numPointBC = 0; ierr = GridAddPointBC(grid, x, y, z, field, f, reduce); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridAddPointBC" /*@C GridAddPointBC This function adds a boundary condition to use for the problem at a point. Collective on Grid Input Parameters: + grid - The grid . x,y,z - The point at which conditions are applied . field - The field to which the boundary condition is applied . f - The function which defines the boundary condition - reduce - The flag for explicit reduction of the system Level: intermediate .keywords active field .seealso GridSetBC, GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction @*/ int GridAddPointBC(Grid grid, double x, double y, double z, int field, PointFunction f, PetscTruth reduce) { GridBC *tempBC; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); GridValidField(grid, field); while (grid->numPointBC >= grid->maxPointBC) { ierr = PetscMalloc(grid->maxPointBC*2 * sizeof(GridBC), &tempBC); CHKERRQ(ierr); PetscLogObjectMemory(grid, grid->maxPointBC * sizeof(GridBC)); ierr = PetscMemcpy(tempBC, grid->pointBC, grid->maxPointBC * sizeof(GridBC)); CHKERRQ(ierr); ierr = PetscFree(grid->pointBC); CHKERRQ(ierr); grid->pointBC = tempBC; grid->maxPointBC *= 2; } if (GridGetNearestBdNode(grid, field, x, y, z, &grid->pointBC[grid->numPointBC].node)) { SETERRQ3(PETSC_ERR_ARG_WRONG, "Invalid point {%g,%g,%g} specified for boundary condition", x, y, z); } grid->pointBC[grid->numPointBC].point[0] = x; grid->pointBC[grid->numPointBC].point[1] = y; grid->pointBC[grid->numPointBC].point[2] = z; grid->pointBC[grid->numPointBC].field = field; grid->pointBC[grid->numPointBC].func = f; grid->pointBC[grid->numPointBC].reduce = reduce; grid->pointBC[grid->numPointBC].boundary = -1; grid->numPointBC++; /* Check whether to reduce system */ if (reduce == PETSC_TRUE) grid->reduceSystem = PETSC_TRUE; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetBCMultiplier" /*@ GridSetBCMultiplier - This sets the scalar multiplier used for reduction components on the rhs. Collective on Grid Input Parameters: + grid - The grid - alpha - The scalar multiplier Note: For example, this should be -1 in a nonlinear iteration. The default is 1. Level: developer .keywords: grid, reduction, boundary conditions .seealso: GridGetBCMultiplier(), GridSetBC(), GridAddBC() @*/ int GridSetBCMultiplier(Grid grid, PetscScalar alpha) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); grid->reduceAlpha = alpha; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetBCMultiplier" /*@ GridGetBCMultiplier - This gets the scalar multiplier used for reduction components on the rhs. Not collective Input Parameter: . grid - The grid Output Parameter: . alpha - The scalar multiplier Level: developer .keywords: grid, reduction, boundary conditions .seealso: GridSetBCMultiplier(), GridSetBC(), GridAddBC() @*/ int GridGetBCMultiplier(Grid grid, PetscScalar *alpha) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidScalarPointer(alpha); *alpha = grid->reduceAlpha; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetBCContext" /*@ GridSetBCContext - This sets the optional user context passed to all routines which assemble boundary reduction information. Must be called before GridSetUp(). Collective on Grid Input Parameters: + grid - The grid - ctx - The context Level: intermediate .keywords: grid, reduction, boundary conditions .seealso: GridGetBCContext(), GridSetBC(), GridAddBC() @*/ int GridSetBCContext(Grid grid, void *ctx) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); grid->reduceContext = ctx; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetBCContext" /*@ GridGetBCContext - This gets the optional user context passed to all routines which assemble boundary reduction information. Not collective Input Parameter: . grid - The grid Output parameter: . ctx - The context Level: intermediate .keywords: grid, reduction, boundary conditions .seealso: GridSetBCContext(), GridSetBC(), GridAddBC() @*/ int GridGetBCContext(Grid grid, void **ctx) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(ctx); *ctx = grid->reduceContext; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetBCValuesType" /*@ GridSetBCValuesType - This determines which boundary values are used to reduce the system. It is intended to allow time dependent boundary conditions to be used, and also supports the difference of two sets of values. Collective on Grid Input Parameter: . grid - The grid Level: intermediate .keywords: grid, reduction, boundary conditions .seealso: GridSetBC(), GridAddBC() @*/ int GridSetBCValuesType(Grid grid, BCValuesType type) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); if (grid->reduceSystem == PETSC_FALSE) PetscFunctionReturn(0); switch(type) { case BC_VALUES: grid->bdReduceVecCur = grid->bdReduceVec; break; case BC_VALUES_OLD: if (grid->bdReduceVecOld == PETSC_NULL) { SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Old boundary values not stored"); } grid->bdReduceVecCur = grid->bdReduceVecOld; break; case BC_VALUES_DIFF: if (grid->bdReduceVecDiff == PETSC_NULL) { SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Difference of boundary values not stored"); } grid->bdReduceVecCur = grid->bdReduceVecDiff; break; default: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid type %d for boundary value calculation", type); } PetscFunctionReturn(0); } /*---------------------------------------------- Calculation Functions ----------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridCalcPointBCNodes" /*@C GridCalcPointBCNodes This function recalculates the nodes used for point boundary conditions. Collective on Grid Input Parameter: . grid - The grid Notes: This function is called by GridReform() after the mesh is recalculated. Level: advanced .keywords grid, point BC, node .seealso GridSetBC, GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction @*/ int GridCalcPointBCNodes(Grid grid) { double x, y, z; int bc; PetscFunctionBegin; for(bc = 0; bc < grid->numPointBC; bc++) { x = grid->pointBC[bc].point[0]; y = grid->pointBC[bc].point[1]; z = grid->pointBC[bc].point[2]; if (GridGetNearestBdNode(grid, grid->pointBC[bc].field, x, y, z, &grid->pointBC[bc].node)) { SETERRQ3(PETSC_ERR_ARG_WRONG, "Invalid point {%g,%g,%g} specified for boundary condition", x, y, z); } } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSaveBCValues_Private" int GridSaveBCValues_Private(Grid grid, VarOrdering reduceOrder, Vec reduceVec) { PetscScalar *array, *arrayOld; int ierr; PetscFunctionBegin; /* Create storage for reduction of Rhs */ if (grid->bdReduceVecOld == PETSC_NULL) { ierr = GVecDuplicate(reduceVec, &grid->bdReduceVecOld); CHKERRQ(ierr); } else if (((Vec_MPI *) grid->bdReduceVecOld->data)->nghost != ((Vec_MPI *) grid->bdReduceVec->data)->nghost) { ierr = GVecDestroy(grid->bdReduceVecOld); CHKERRQ(ierr); ierr = GVecDuplicate(reduceVec, &grid->bdReduceVecOld); CHKERRQ(ierr); } /* ierr = VecCopy(grid->bdReduceVec, grid->bdReduceVecOld); CHKERRQ(ierr); */ ierr = VecGetArray(reduceVec, &array); CHKERRQ(ierr); ierr = VecGetArray(grid->bdReduceVecOld, &arrayOld); CHKERRQ(ierr); ierr = PetscMemcpy(arrayOld, array, reduceOrder->numOverlapVars * sizeof(PetscScalar)); CHKERRQ(ierr); ierr = VecRestoreArray(reduceVec, &array); CHKERRQ(ierr); ierr = VecRestoreArray(grid->bdReduceVecOld, &arrayOld); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridCalcGridBCValues_Private" int GridCalcGridBCValues_Private(Grid grid, VarOrdering reduceOrder, Vec reduceVec, void *ctx) { GridBC *gBC = grid->bc; VarOrdering bcOrder; int bc; int ierr; PetscFunctionBegin; /* Evaluate the vector of boundary values -- If order->localStart[field] is NULL, this means the field is not present in the ordering. This is a better check than seeing if the field is active, since we might want to pass in an order on that field to make boundary values for an inactive field. */ for(bc = 0; bc < grid->numBC; bc++) { if (gBC[bc].reduce != PETSC_TRUE) continue; if (reduceOrder->localStart[gBC[bc].field] == PETSC_NULL) continue; ierr = VarOrderingCreateSubset(reduceOrder, 1, &gBC[bc].field, PETSC_FALSE, &bcOrder); CHKERRQ(ierr); ierr = (*grid->ops->gvecevaluatefunctionboundary)(grid, reduceVec, gBC[bc].boundary, bcOrder, gBC[bc].func, 1.0, ctx);CHKERRQ(ierr); ierr = VarOrderingDestroy(bcOrder); CHKERRQ(ierr); #ifdef PETSC_USE_BOPT_g ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridCalcPointBCValues_Private" int GridCalcPointBCValues_Private(Grid grid, VarOrdering reduceOrder, Vec reduceVec, void *ctx) { GridBC *pBC = grid->pointBC; int **localStart = reduceOrder->localStart; int *offsets = reduceOrder->offsets; int *localOffsets = reduceOrder->localOffsets; FieldClassMap map; VarOrdering bcOrder; PetscScalar *array; double x, y, z; int numNodes, firstVar, rank; int bc, field, numComp, node, nclass, row; int ierr; PetscFunctionBegin; ierr = MPI_Comm_rank(grid->comm, &rank); CHKERRQ(ierr); ierr = VarOrderingGetClassMap(reduceOrder, &map); CHKERRQ(ierr); numNodes = map->numNodes; firstVar = reduceOrder->firstVar[rank]; /* Evaluate the vector of boundary values -- If order->localStart[field] is NULL, this means the field is not present in the ordering. This is a better check than seeing if the field is active, since we might want to pass in an order on that field to make boundary values for an inactive field. */ ierr = VecGetArray(reduceVec, &array); CHKERRQ(ierr); for(bc = 0; bc < grid->numPointBC; bc++) { /* BC is not reduced out of the system */ if (pBC[bc].reduce != PETSC_TRUE) continue; /* Field is not present in the ordering */ if (reduceOrder->localStart[pBC[bc].field] == PETSC_NULL) continue; field = pBC[bc].field; numComp = grid->fields[field].numComp; ierr = VarOrderingCreateSubset(reduceOrder, 1, &field, PETSC_FALSE, &bcOrder); CHKERRQ(ierr); /* Point is in another domain */ if (pBC[bc].node < 0) continue; node = pBC[bc].node; nclass = map->classes[node]; if (node >= numNodes) { row = localOffsets[node-numNodes]; } else { row = offsets[node] - firstVar + localStart[field][nclass]; } x = pBC[bc].point[0]; y = pBC[bc].point[1]; z = pBC[bc].point[2]; ierr = (*pBC[bc].func)(1, numComp, &x, &y, &z, &array[row], ctx); CHKERRQ(ierr); ierr = VarOrderingDestroy(bcOrder); CHKERRQ(ierr); #ifdef PETSC_USE_BOPT_g ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif } ierr = VecRestoreArray(reduceVec, &array); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridCalcBCValues_Private" int GridCalcBCValues_Private(Grid grid, VarOrdering reduceOrder, Vec reduceVec, PetscTruth save, void *ctx) { PetscScalar *array; int numGhostVars; int ierr; PetscFunctionBegin; numGhostVars = reduceOrder->numOverlapVars - reduceOrder->numLocVars; if (((Vec_MPI *) reduceVec->data)->nghost != numGhostVars) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid reduce vector size %d should be %d", ((Vec_MPI *) reduceVec->data)->nghost, numGhostVars); } if (save == PETSC_TRUE) { ierr = GridSaveBCValues_Private(grid, reduceOrder, reduceVec); CHKERRQ(ierr); } /* Initialize vector */ /* ierr = VecSet(&zero, reduceVec); CHKERRQ(ierr); */ ierr = VecGetArray(reduceVec, &array); CHKERRQ(ierr); ierr = PetscMemzero(array, reduceOrder->numOverlapVars * sizeof(PetscScalar)); CHKERRQ(ierr); ierr = VecRestoreArray(reduceVec, &array); CHKERRQ(ierr); ierr = GridCalcGridBCValues_Private(grid, reduceOrder, reduceVec, ctx); CHKERRQ(ierr); ierr = GridCalcPointBCValues_Private(grid, reduceOrder, reduceVec, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridCalcBCValues" /*@ GridCalcBCValues - This function calculates the boundary values. It is normally called once a timestep when using time dependent boundary conditions. Collective on Grid Input Parameters: + grid - The grid . save - A flag used to store old values, usually for timestepping - ctx - The context Level: advanced .keywords: grid, reduction, boundary conditions .seealso: GridSetBCContext(), GridSetBC(), GridAddBC() @*/ int GridCalcBCValues(Grid grid, PetscTruth save, void *ctx) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); if (grid->reduceSystem == PETSC_TRUE) { if (ctx != PETSC_NULL) PetscValidPointer(ctx); ierr = GridCalcBCValues_Private(grid, grid->reduceOrder, grid->bdReduceVec, save, ctx); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridCalcBCValuesDifference" /*@ GridCalcBCValuesDifference - This function calculates the difference of the last two sets of boundary values and puts it in an internal vector. This is is normally used to implement the Rhs time derivative in a GTS. Collective on Grid Input Parameter: . grid - The grid Level: advanced .keywords: grid, reduction, boundary conditions .seealso: GridSetBCContext(), GridSetBC(), GridAddBC() @*/ int GridCalcBCValuesDifference(Grid grid) { PetscScalar *array, *arrayOld, *arrayDiff; int numGhostVars; register int i, n; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); if (grid->reduceSystem == PETSC_TRUE) { numGhostVars = grid->reduceOrder->numOverlapVars - grid->reduceOrder->numLocVars; if (((Vec_MPI *) grid->bdReduceVec->data)->nghost != numGhostVars) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduce vector size"); } if (grid->bdReduceVecOld == PETSC_NULL) { SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "No previous boundary values"); } /* Create storage for reduction of Rhs */ if (grid->bdReduceVecDiff == PETSC_NULL) { ierr = GVecDuplicate(grid->bdReduceVec, &grid->bdReduceVecDiff); CHKERRQ(ierr); } else if (((Vec_MPI *) grid->bdReduceVecDiff->data)->nghost != ((Vec_MPI *) grid->bdReduceVec->data)->nghost) { ierr = GVecDestroy(grid->bdReduceVecDiff); CHKERRQ(ierr); ierr = GVecDuplicate(grid->bdReduceVec, &grid->bdReduceVecDiff); CHKERRQ(ierr); } /* ierr = VecWAXPY(&minusOne, grid->bdReduceVecOld, grid->bdReduceVec, grid->bdReduceVecDiff); CHKERRQ(ierr); */ ierr = VecGetArray(grid->bdReduceVec, &array); CHKERRQ(ierr); ierr = VecGetArray(grid->bdReduceVecOld, &arrayOld); CHKERRQ(ierr); ierr = VecGetArray(grid->bdReduceVecDiff, &arrayDiff); CHKERRQ(ierr); n = grid->reduceOrder->numOverlapVars; PetscLogFlops(n); for(i = 0; i < n; i++) arrayDiff[i] = array[i] - arrayOld[i]; ierr = VecRestoreArray(grid->bdReduceVec, &array); CHKERRQ(ierr); ierr = VecRestoreArray(grid->bdReduceVecOld, &arrayOld); CHKERRQ(ierr); ierr = VecRestoreArray(grid->bdReduceVecDiff, &arrayDiff); CHKERRQ(ierr); } #ifdef PETSC_USE_BOPT_g ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif PetscFunctionReturn(0); } /*----------------------------------------------- Reduction Functions -----------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridSetReduceSystem" /*@C GridSetReduceSystem - This function determines whether unknowns associated with boundary conditions are eliminated from the system. Collective on Grid Input Parameters: + grid - The grid - reduce - The flag for explicit reduction of the system Level: intermediate .keywords grid, boundary condition, reduce .seealso GridGetReduceSystem(), GridSetBC(), GridAddBC(), GridSetPointBC(), GridAddPointBC() @*/ int GridSetReduceSystem(Grid grid, PetscTruth reduce) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); grid->reduceSystem = reduce; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetReduceSystem" /*@C GridGetReduceSystem - This function reveals whether unknowns associated with boundary conditions are eliminated from the system. Not collective Input Parameter: . grid - The grid Output Parameter: . reduce - The flag for explicit reduction of the system Level: intermediate .keywords grid, boundary condition, reduce .seealso GridSetReduceSystem(), GridSetBC(), GridAddBC(), GridSetPointBC(), GridAddPointBC() @*/ int GridGetReduceSystem(Grid grid, PetscTruth *reduce) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(reduce); *reduce = grid->reduceSystem; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetReduceElement" /*@C GridSetReduceElement - This function determines whether element matrices and vectors are reduced on the fly, or if boundary operators are stored and applied. Collective on Grid Input Parameters: + grid - The grid - reduce - The flag for explicit reduction of the system Level: intermediate .keywords grid, boundary condition, reduce, element .seealso GridGetReduceElement(), GridSetBC(), GridAddBC(), GridSetPointBC(), GridAddPointBC() @*/ int GridSetReduceElement(Grid grid, PetscTruth reduce) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); grid->reduceElement = reduce; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridGetReduceElement" /*@C GridGetReduceElement - This function indicates whether element matrices and vectors are reduced on the fly, or if boundary operators are stored and applied. Not collective Input Parameter: . grid - The grid Output Parameter: . reduce - The flag for explicit reduction of the system Level: intermediate .keywords grid, boundary condition, reduce, element .seealso GridSetReduceElement(), GridSetBC(), GridAddBC(), GridSetPointBC(), GridAddPointBC() @*/ int GridGetReduceElement(Grid grid, PetscTruth *reduce) { PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(reduce); *reduce = grid->reduceElement; PetscFunctionReturn(0); } /*---------------------------------------------- Application Functions ----------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridResetConstrainedMultiply_Private" /* GridResetConstrainedMultiply_Private - This function resets the mulplication routine for constrained matrices Input Parameters: + grid - The Grid - A - The GMat Level: developer .keywords Grid, GMat, reset, constrained, multiply .seealso GridEvaluateRhs */ int GridResetConstrainedMultiply_Private(Grid grid, GMat A) { PetscTruth isConstrained, explicitConstraints; void (*oldMult)(void); int ierr; PetscFunctionBegin; ierr = GridIsConstrained(grid, &isConstrained); CHKERRQ(ierr); ierr = GridGetExplicitConstraints(grid, &explicitConstraints); CHKERRQ(ierr); if (isConstrained == PETSC_TRUE) { if (explicitConstraints == PETSC_FALSE) { ierr = MatShellGetOperation(A, MATOP_MULT, &oldMult); CHKERRQ(ierr); if (oldMult != (void (*)(void)) GMatMatMultConstrained) { ierr = MatShellSetOperation(A, MATOP_MULT_CONSTRAINED, oldMult); CHKERRQ(ierr); } ierr = MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) GMatMatMultConstrained); CHKERRQ(ierr); ierr = MatShellGetOperation(A, MATOP_MULT_TRANSPOSE, &oldMult); CHKERRQ(ierr); if (oldMult != (void (*)(void)) GMatMatMultTransposeConstrained) { ierr = MatShellSetOperation(A, MATOP_MULT_TRANSPOSE_CONSTRAINED, oldMult); CHKERRQ(ierr); } ierr = MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, (void (*)(void)) GMatMatMultTransposeConstrained); CHKERRQ(ierr); } else { ierr = MatShellGetOperation(A, MATOP_MULT_CONSTRAINED, &oldMult); CHKERRQ(ierr); if (oldMult != PETSC_NULL) { ierr = MatShellSetOperation(A, MATOP_MULT, oldMult); CHKERRQ(ierr); } ierr = MatShellGetOperation(A, MATOP_MULT_TRANSPOSE_CONSTRAINED, &oldMult); CHKERRQ(ierr); if (oldMult != PETSC_NULL) { ierr = MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, oldMult); CHKERRQ(ierr); } } } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetBoundary" /*@C GridSetBoundary This function sets Dirchlet boundary conditions on the linear problem arising from the underlying grid. Collective on Grid Input Parameters: + bd - The marker for the boundary to apply conditions along . field - The field to which the conditions apply . diag - The scalar to be placed on the diagonal . f - The function which defines the boundary condition - ctx - The user-supplied context Output Parameters: + A - The system matrix - b - The Rhs vector Level: intermediate .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetBoundary(int bd, int field, PetscScalar diag, PointFunction f, GMat A, GVec b, void *ctx) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(A, MAT_COOKIE); ierr = GMatGetGrid(A, &grid); CHKERRQ(ierr); ierr = GridSetBoundaryRectangular(bd, field, diag, f, grid->order, A, b, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetBoundaryRectangular" /*@C GridSetBoundaryRectangular This function sets Dirchlet boundary conditions on the linear problem arising from the underlying grid, and the default variable ordering can be overridden. Collective on Grid Input Parameters: + bd - The marker for the boundary to apply conditions along . field - The field to which the conditions apply . diag - The scalar to be placed on the diagonal . f - The function which defines the boundary condition . order - The test variable ordering - ctx - The user-supplied context Output Parameters: + A - The system matrix - b - The Rhs vector Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetBoundaryRectangular(int bd, int field, PetscScalar diag, PointFunction f, VarOrdering order, GMat A, GVec b, void *ctx) { Grid grid, grid2; Mesh mesh; int comp; /* The number of field components */ int size; /* The number of nodes in the boundary */ int *localStart; /* The offset of this field on a node of a given class */ int node; /* The canonical node number of the current boundary node */ int nclass; /* The class of the current boundary node */ double *x, *y, *z; /* Coordinates of the boundary nodes */ int vars; /* The number of variables affected (var/node * size) */ int *offsets; /* The canonical variable number for the first variable on each node */ int *rows; /* Rows corresponding to boundary nodes */ PetscScalar *values; /* Boundary values */ PetscScalar elem = diag; IS is; int rank; int i, j, count; #ifdef PETSC_USE_BOPT_g PetscTruth opt; #endif int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(A, MAT_COOKIE); PetscValidHeaderSpecific(b, VEC_COOKIE); PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE); ierr = GMatGetGrid(A, &grid); CHKERRQ(ierr); ierr = GVecGetGrid(b, &grid2); CHKERRQ(ierr); if (grid != grid2) SETERRQ(PETSC_ERR_ARG_INCOMP, "Matrix and vector have different underlying grids"); GridValidField(grid, field); ierr = MPI_Comm_rank(grid->comm, &rank); CHKERRQ(ierr); mesh = grid->mesh; comp = grid->fields[field].disc->comp; offsets = order->offsets; localStart = order->localStart[field]; /* Support for constrained problems */ ierr = VecGetSize(b, &size); CHKERRQ(ierr); if (grid->isConstrained) { if (size != grid->constraintOrder->numVars) SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid vector size"); offsets = grid->constraintOrder->offsets; } else { if (size != grid->order->numVars) SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid vector size"); } /* Allocate memory */ ierr = GridGetBoundarySize(grid, bd, field, &size); CHKERRQ(ierr); if (size == 0) { #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscSynchronizedFlush(grid->comm); PetscSynchronizedFlush(grid->comm); } #endif ierr = VecAssemblyBegin(b); CHKERRQ(ierr); ierr = VecAssemblyEnd(b); CHKERRQ(ierr); ierr = ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is); CHKERRQ(ierr); ierr = MatZeroRows(A, is, &elem); CHKERRQ(ierr); ierr = ISDestroy(is); CHKERRQ(ierr); PetscFunctionReturn(0); } vars = size*comp; ierr = PetscMalloc(size * sizeof(double), &x); CHKERRQ(ierr); ierr = PetscMalloc(size * sizeof(double), &y); CHKERRQ(ierr); ierr = PetscMalloc(size * sizeof(double), &z); CHKERRQ(ierr); ierr = PetscMalloc(vars * sizeof(PetscScalar), &values); CHKERRQ(ierr); ierr = PetscMalloc(vars * sizeof(int), &rows); CHKERRQ(ierr); /* Loop over boundary nodes */ ierr = GridGetBoundaryStart(grid, bd, field, PETSC_FALSE, &node, &nclass); CHKERRQ(ierr); for(i = 0, count = 0; node >= 0; i++) { for(j = 0; j < comp; j++, count++) { rows[count] = offsets[node] + j + localStart[nclass]; #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscSynchronizedPrintf(grid->comm, "[%d]bd %d field: %d node: %d row: %d class: %d\n", rank, bd, field, node, rows[count], nclass); } #endif } ierr = MeshGetNodeCoords(mesh, node, &x[i], &y[i], &z[i]); CHKERRQ(ierr); ierr = GridGetBoundaryNext(grid, bd, field, PETSC_FALSE, &node, &nclass); CHKERRQ(ierr); } #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscSynchronizedFlush(grid->comm); } #endif /* Get boundary values */ ierr = (*f)(size, comp, x, y, z, values, ctx); CHKERRQ(ierr); /* Put values in Rhs */ #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscPrintf(grid->comm, "Setting boundary values in rhs bd %d field %d\n", bd, field); for(i = 0; i < vars; i++) PetscSynchronizedPrintf(grid->comm, " row: %d val: %g\n", rows[i], PetscRealPart(values[i])); PetscSynchronizedFlush(grid->comm); } ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif ierr = VecSetValues(b, vars, rows, values, INSERT_VALUES); CHKERRQ(ierr); ierr = VecAssemblyBegin(b); CHKERRQ(ierr); ierr = VecAssemblyEnd(b); CHKERRQ(ierr); /* Set rows of A to the identity */ ierr = ISCreateGeneral(PETSC_COMM_SELF, vars, rows, &is); CHKERRQ(ierr); ierr = MatZeroRows(A, is, &elem); CHKERRQ(ierr); ierr = ISDestroy(is); CHKERRQ(ierr); ierr = GridResetConstrainedMultiply_Private(grid, A); CHKERRQ(ierr); ierr = PetscFree(x); CHKERRQ(ierr); ierr = PetscFree(y); CHKERRQ(ierr); ierr = PetscFree(z); CHKERRQ(ierr); ierr = PetscFree(values); CHKERRQ(ierr); ierr = PetscFree(rows); CHKERRQ(ierr); PetscFunctionReturn(0); } /*------------------------------------------------- Matrix Functions ------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridSetMatBoundary" /*@C GridSetMatBoundary This function sets Dirchlet boundary conditions on the linear system matrix arising from the underlying grid. Collective on GMat Input Parameters: + bd - The marker for the boundary to apply conditions along . field - The field to which the conditions apply . diag - The scalar to be placed on the diagonal - ctx - The user-supplied context Output Parameter: . A - The system matrix Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetMatBoundary(int bd, int field, PetscScalar diag, GMat A, void *ctx) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(A, MAT_COOKIE); ierr = GMatGetGrid(A, &grid); CHKERRQ(ierr); ierr = GridSetMatBoundaryRectangular(1, &bd, &field, diag, grid->order, A, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetMatBoundaryRectangular" /*@C GridSetMatBoundaryRectangular This function sets Dirchlet boundary conditions on the linear system matrix arising from the underlying grid, and the default variable ordering can be overridden. Collective on GMat Input Parameters: + num - The number of boundary conditions . bd - The markers for each boundary to apply conditions along . field - The fields to which the conditions apply . diag - The scalar to be placed on the diagonal . order - The test variable ordering - ctx - The user-supplied context Output Parameter: . A - The system matrix Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetMatBoundaryRectangular(int num, int *bd, int *field, PetscScalar diag, VarOrdering order, GMat A, void *ctx) { Grid grid; int comp; /* The number of field components */ int size; /* The number of nodes in the boundary */ int totSize; /* The number of nodes in all boundaries */ int *localStart; /* The offset of this field on a node of a given class */ int node; /* The canonical node number of the current boundary node */ int nclass; /* The class of the current boundary node */ int vars; /* The number of variables affected (var/node * size) */ int *offsets; /* The canonical variable number for the first variable on each node */ int *rows; /* Rows corresponding to boundary nodes */ PetscScalar elem = diag; IS is; int rank; int b, i, j, count; #ifdef PETSC_USE_BOPT_g PetscTruth opt; #endif int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(A, MAT_COOKIE); PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE); ierr = GMatGetGrid(A, &grid); CHKERRQ(ierr); offsets = order->offsets; ierr = MPI_Comm_rank(grid->comm, &rank); CHKERRQ(ierr); /* Allocate memory */ for(b = 0, totSize = 0, vars = 0; b < num; b++) { GridValidField(grid, field[b]); ierr = GridGetBoundarySize(grid, bd[b], field[b], &size); CHKERRQ(ierr); totSize += size; vars += size*grid->fields[field[b]].disc->comp; } if (totSize == 0) { #ifdef PETSC_USE_BOPT_g PetscSynchronizedFlush(grid->comm); #endif ierr = ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is); CHKERRQ(ierr); ierr = MatZeroRows(A, is, &elem); CHKERRQ(ierr); ierr = ISDestroy(is); CHKERRQ(ierr); PetscFunctionReturn(0); } ierr = PetscMalloc(vars * sizeof(int), &rows); CHKERRQ(ierr); /* Loop over boundaries */ for(b = 0, count = 0; b < num; b++) { comp = grid->fields[field[b]].disc->comp; localStart = order->localStart[field[b]]; /* Loop over boundary nodes */ ierr = GridGetBoundaryStart(grid, bd[b], field[b], PETSC_FALSE, &node, &nclass); CHKERRQ(ierr); for(i = 0; node >= 0; i++) { for(j = 0; j < comp; j++, count++) { rows[count] = offsets[node] + j + localStart[nclass]; #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscSynchronizedPrintf(grid->comm, "[%d]bd %d field: %d node: %d row: %d class: %d\n", rank, bd[b], field[b], node, rows[count], nclass); } #endif } ierr = GridGetBoundaryNext(grid, bd[b], field[b], PETSC_FALSE, &node, &nclass); CHKERRQ(ierr); } } #ifdef PETSC_USE_BOPT_g PetscSynchronizedFlush(grid->comm); if (count != vars) SETERRQ2(PETSC_ERR_PLIB, "Boundary size %d should be %d", count, vars); ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif /* Set rows of A to the identity */ ierr = ISCreateGeneral(PETSC_COMM_SELF, vars, rows, &is); CHKERRQ(ierr); ierr = MatZeroRows(A, is, &elem); CHKERRQ(ierr); ierr = ISDestroy(is); CHKERRQ(ierr); ierr = GridResetConstrainedMultiply_Private(grid, A); CHKERRQ(ierr); ierr = PetscFree(rows); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetMatPointBoundary" /*@C GridSetMatPointBoundary This function sets Dirchlet boundary conditions on the linear system matrix arising from the underlying grid. Collective on GMat Input Parameters: + node - The constrained node . field - The field to which the conditions apply . diag - The scalar to be placed on the diagonal - ctx - The user-supplied context Output Parameter: . A - The system matrix Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetMatPointBoundary(int node, int field, PetscScalar diag, GMat A, void *ctx) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(A, MAT_COOKIE); ierr = GMatGetGrid(A, &grid); CHKERRQ(ierr); ierr = GridSetMatPointBoundaryRectangular(node, field, diag, grid->order, A, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetMatPointBoundaryRectangular" /*@C GridSetMatPointBoundaryRectangular This function sets Dirchlet boundary conditions on the linear system matrix arising from the underlying grid, and the default variable ordering can be overridden. Collective on GMat Input Parameters: + node - The constrained node . field - The field to which the conditions apply . diag - The scalar to be placed on the diagonal . order - The test variable ordering - ctx - The user-supplied context Output Parameter: . A - The system matrix Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetMatPointBoundaryRectangular(int node, int field, PetscScalar diag, VarOrdering order, GMat A, void *ctx) { Grid grid; int comp; /* The number of field components */ int *localStart; /* The offset of this field on a node of a given class */ int nclass; /* The class of the current boundary node */ int *offsets; /* The canonical variable number for the first variable on each node */ int *rows; /* Rows corresponding to boundary nodes */ PetscScalar elem = diag; IS is; int rank; int j; #ifdef PETSC_USE_BOPT_g PetscTruth opt; #endif int ierr; PetscFunctionBegin; if (node < 0) { ierr = ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is); CHKERRQ(ierr); ierr = MatZeroRows(A, is, &elem); CHKERRQ(ierr); ierr = ISDestroy(is); CHKERRQ(ierr); PetscFunctionReturn(0); } PetscValidHeaderSpecific(A, MAT_COOKIE); PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE); ierr = GMatGetGrid(A, &grid); CHKERRQ(ierr); GridValidField(grid, field); ierr = MPI_Comm_rank(grid->comm, &rank); CHKERRQ(ierr); comp = grid->fields[field].disc->comp; offsets = order->offsets; localStart = order->localStart[field]; /* Allocate memory */ ierr = PetscMalloc(comp * sizeof(int), &rows); CHKERRQ(ierr); ierr = GridGetNodeClass(grid, node, &nclass); CHKERRQ(ierr); for(j = 0; j < comp; j++) { rows[j] = offsets[node] + j + localStart[nclass]; #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscPrintf(PETSC_COMM_SELF, "[%d]field: %d node: %d row: %d class: %d\n", rank, field, node, rows[j], nclass); } #endif } #ifdef PETSC_USE_BOPT_g ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif /* Set rows of A to the identity */ ierr = ISCreateGeneral(PETSC_COMM_SELF, comp, rows, &is); CHKERRQ(ierr); ierr = MatZeroRows(A, is, &elem); CHKERRQ(ierr); ierr = ISDestroy(is); CHKERRQ(ierr); ierr = GridResetConstrainedMultiply_Private(grid, A); CHKERRQ(ierr); ierr = PetscFree(rows); CHKERRQ(ierr); PetscFunctionReturn(0); } /*------------------------------------------------- Vector Functions ------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridSetVecBoundary" /*@C GridSetVecBoundary This function sets Dirchlet boundary conditions on the linear Rhs arising from the underlying grid. Collective on GVec Input Parameters: + bd - The marker for the boundary to apply conditions along . field - The field to which the conditions apply . f - The function which defines the boundary condition - ctx - The user-supplied context Output Parameter: . b - The Rhs vector Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetVecBoundary(int bd, int field, PointFunction f, GVec b, void *ctx) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(b, VEC_COOKIE); ierr = GVecGetGrid(b, &grid); CHKERRQ(ierr); ierr = GridSetVecBoundaryRectangular(1, &bd, &field, &f, grid->order, b, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetVecBoundaryRectangular" /*@C GridSetVecBoundaryRectangular This function sets Dirchlet boundary conditions on the linear Rhs arising from the underlying grid, and the default variable ordering can be overridden. Collective on GVec Input Parameters: + num - The number of boundary conditions . bd - The markers for each boundary to apply conditions along . field - The fields to which the conditions apply . f - The functions which define the boundary conditions . order - The test variable ordering - ctx - The user-supplied context Output Parameter: . b - The Rhs vector Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetVecBoundaryRectangular(int num, int *bd, int *field, PointFunction *f, VarOrdering order, GVec b, void *ctx) { Grid grid; Mesh mesh; int comp; /* The number of field components */ int *sizes; /* The number of nodes in each boundary */ int totSize; /* The number of nodes in all boundaries */ int maxSize; /* The maximum number of nodes in any boundary */ int *localStart; /* The offset of this field on a node of a given class */ int node; /* The canonical node number of the current boundary node */ int nclass; /* The class of the current boundary node */ double *x, *y, *z; /* Coordinates of the boundary nodes */ int vars; /* The number of variables affected (var/node * size) */ int *offsets; /* The canonical variable number for the first variable on each node */ int *rows; /* Rows corresponding to boundary nodes */ PetscScalar *values; /* Boundary values */ int size, rank; int c, i, j, count, off; #ifdef PETSC_USE_BOPT_g PetscTruth opt; #endif int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(b, VEC_COOKIE); PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE); ierr = GVecGetGrid(b, &grid); CHKERRQ(ierr); mesh = grid->mesh; offsets = order->offsets; ierr = MPI_Comm_rank(grid->comm, &rank); CHKERRQ(ierr); /* Support for constrained problems */ ierr = VecGetSize(b, &size); CHKERRQ(ierr); if (grid->isConstrained) { if (size != grid->constraintOrder->numVars) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->constraintOrder->numVars); } offsets = grid->constraintOrder->offsets; } else { if (size != grid->order->numVars) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->order->numVars); } } /* Allocate memory */ ierr = PetscMalloc(num * sizeof(int), &sizes); CHKERRQ(ierr); for(c = 0, totSize = 0, maxSize = 0, vars = 0; c < num; c++) { GridValidField(grid, field[c]); ierr = GridGetBoundarySize(grid, bd[c], field[c], &sizes[c]); CHKERRQ(ierr); totSize += sizes[c]; maxSize = PetscMax(maxSize, sizes[c]); vars += sizes[c]*grid->fields[field[c]].disc->comp; } if (totSize == 0) { #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscSynchronizedFlush(grid->comm); PetscSynchronizedFlush(grid->comm); } #endif ierr = VecAssemblyBegin(b); CHKERRQ(ierr); ierr = VecAssemblyEnd(b); CHKERRQ(ierr); PetscFunctionReturn(0); } ierr = PetscMalloc(maxSize * sizeof(double), &x); CHKERRQ(ierr); ierr = PetscMalloc(maxSize * sizeof(double), &y); CHKERRQ(ierr); ierr = PetscMalloc(maxSize * sizeof(double), &z); CHKERRQ(ierr); ierr = PetscMalloc(vars * sizeof(PetscScalar), &values); CHKERRQ(ierr); ierr = PetscMalloc(vars * sizeof(int), &rows); CHKERRQ(ierr); /* Loop over boundaries */ for(c = 0, count = 0, off = 0; c < num; c++, off = count) { if (sizes[c] == 0) { #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscSynchronizedFlush(grid->comm); PetscSynchronizedFlush(grid->comm); } #endif continue; } comp = grid->fields[field[c]].disc->comp; localStart = order->localStart[field[c]]; /* Loop over boundary nodes */ ierr = GridGetBoundaryStart(grid, bd[c], field[c], PETSC_FALSE, &node, &nclass); CHKERRQ(ierr); for(i = 0; node >= 0; i++) { for(j = 0; j < comp; j++, count++) { rows[count] = offsets[node] + j + localStart[nclass]; #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscSynchronizedPrintf(grid->comm, "[%d]bd %d field: %d node: %d row: %d class: %d\n", rank, bd[c], field[c], node, rows[count], nclass); } #endif } ierr = MeshGetNodeCoords(mesh, node, &x[i], &y[i], &z[i]); CHKERRQ(ierr); ierr = GridGetBoundaryNext(grid, bd[c], field[c], PETSC_FALSE, &node, &nclass); CHKERRQ(ierr); } #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscSynchronizedFlush(grid->comm); } #endif /* Get boundary values */ ierr = (*(f[c]))(sizes[c], comp, x, y, z, &values[off], ctx); CHKERRQ(ierr); #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscPrintf(grid->comm, "Setting boundary values in rhs bd %d field %d\n", bd[c], field[c]); for(i = off; i < count; i++) PetscSynchronizedPrintf(grid->comm, " row: %d val: %g\n", rows[i], PetscRealPart(values[i])); PetscSynchronizedFlush(grid->comm); } ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif } if (count != vars) SETERRQ2(PETSC_ERR_PLIB, "Boundary size %d should be %d", count, vars); /* Put values in Rhs */ ierr = VecSetValues(b, vars, rows, values, INSERT_VALUES); CHKERRQ(ierr); ierr = VecAssemblyBegin(b); CHKERRQ(ierr); ierr = VecAssemblyEnd(b); CHKERRQ(ierr); ierr = PetscFree(sizes); CHKERRQ(ierr); ierr = PetscFree(x); CHKERRQ(ierr); ierr = PetscFree(y); CHKERRQ(ierr); ierr = PetscFree(z); CHKERRQ(ierr); ierr = PetscFree(values); CHKERRQ(ierr); ierr = PetscFree(rows); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetVecPointBoundary" /*@C GridSetVecPointBoundary This function sets Dirchlet boundary conditions on the linear Rhs arising from the underlying grid. Collective on GVec Input Parameters: + node - The constrained node . field - The field to which the conditions apply . f - The function which defines the boundary condition - ctx - The user-supplied context Output Parameter: . b - The Rhs vector Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetVecPointBoundary(int node, int field, PointFunction f, GVec b, void *ctx) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(b, VEC_COOKIE); ierr = GVecGetGrid(b, &grid); CHKERRQ(ierr); ierr = GridSetVecPointBoundaryRectangular(node, field, f, grid->order, b, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetVecPointBoundaryRectangular" /*@C GridSetVecPointBoundaryRectangular This function sets Dirchlet boundary conditions on the linear Rhs arising from the underlying grid, and the default variable ordering can be overridden. Collective on GVec Input Parameters: + node - The constriained node . field - The field to which the conditions apply . f - The function which defines the boundary condition . order - The test variable ordering - ctx - The user-supplied context Output Parameter: . b - The Rhs vector Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetVecPointBoundaryRectangular(int node, int field, PointFunction f, VarOrdering order, GVec b, void *ctx) { Grid grid; Mesh mesh; int comp; /* The number of field components */ int size; /* The number of nodes in the boundary */ int *localStart; /* The offset of this field on a node of a given class */ int nclass; /* The class of the current boundary node */ double x, y, z; /* Coordinates of the boundary nodes */ int *offsets; /* The canonical variable number for the first variable on each node */ int *rows; /* Rows corresponding to boundary nodes */ PetscScalar *values; /* Boundary values */ int rank; int c; #ifdef PETSC_USE_BOPT_g PetscTruth opt; #endif int ierr; PetscFunctionBegin; if (node < 0) { ierr = VecAssemblyBegin(b); CHKERRQ(ierr); ierr = VecAssemblyEnd(b); CHKERRQ(ierr); PetscFunctionReturn(0); } PetscValidHeaderSpecific(b, VEC_COOKIE); PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE); ierr = GVecGetGrid(b, &grid); CHKERRQ(ierr); GridValidField(grid, field); ierr = MPI_Comm_rank(grid->comm, &rank); CHKERRQ(ierr); mesh = grid->mesh; comp = grid->fields[field].disc->comp; offsets = order->offsets; localStart = order->localStart[field]; PetscValidPointer(localStart); /* Support for constrained problems */ ierr = VecGetSize(b, &size); CHKERRQ(ierr); if (grid->isConstrained) { if (size != grid->constraintOrder->numVars) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->constraintOrder->numVars); } offsets = grid->constraintOrder->offsets; } else { if (size != grid->order->numVars) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->order->numVars); } } /* Allocate memory */ size = 1; ierr = PetscMalloc(comp * sizeof(PetscScalar), &values); CHKERRQ(ierr); ierr = PetscMalloc(comp * sizeof(int), &rows); CHKERRQ(ierr); ierr = MeshGetNodeCoords(mesh, node, &x, &y, &z); CHKERRQ(ierr); ierr = GridGetNodeClass(grid, node, &nclass); CHKERRQ(ierr); for(c = 0; c < comp; c++) { rows[c] = offsets[node] + c + localStart[nclass]; #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscPrintf(PETSC_COMM_SELF, "[%d]field: %d node: %d row: %d class: %d\n", rank, field, node, rows[c], nclass); } #endif } /* Get boundary values */ ierr = (*f)(size, comp, &x, &y, &z, values, ctx); CHKERRQ(ierr); /* Put values in Rhs */ #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscPrintf(PETSC_COMM_SELF, "Setting boundary values on rhs node %d field %d\n", node, field); for(c = 0; c < comp; c++) PetscPrintf(PETSC_COMM_SELF, " row: %d val: %g\n", rows[c], PetscRealPart(values[c])); } ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif ierr = VecSetValues(b, comp, rows, values, INSERT_VALUES); CHKERRQ(ierr); ierr = VecAssemblyBegin(b); CHKERRQ(ierr); ierr = VecAssemblyEnd(b); CHKERRQ(ierr); ierr = PetscFree(values); CHKERRQ(ierr); ierr = PetscFree(rows); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetVecBoundaryDifference" /*@C GridSetVecBoundaryDifference This function sets Dirchlet boundary conditions on the linear Rhs arising from the underlying grid, but actually sets it to the difference of the function value and the value in the given vector. This is commonly used in a time dependent, nonlinear problem for which we would like the rhs boundary values to be: U^{n+1}_k - U^{n+1}_{k+1} where n is the time iteration index, and k is the Newton iteration index. This means that the solution will be updated to U^{n+1}_{k+1} if the Jacobian is the identity for that row. This is very useful for time dependent boundary conditions for which the traditional method of letting the rhs value be zero does not work. Collective on GVec Input Parameters: + bd - The marker for the boundary to apply conditions along . u - A grid vector, usually the previous solution . field - The field to which the conditions apply . f - The function which defines the boundary condition - ctx - The user-supplied context Output Parameter: . b - The Rhs vector Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetVecBoundaryDifference(int bd, int field, GVec u, PointFunction f, GVec b, void *ctx) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(b, VEC_COOKIE); ierr = GVecGetGrid(b, &grid); CHKERRQ(ierr); ierr = GridSetVecBoundaryDifferenceRectangular(bd, field, u, f, grid->order, b, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetVecBoundaryDifferenceRectangular" /*@C GridSetVecBoundaryDifferenceRectangular This function sets Dirchlet boundary conditions on the linear Rhs arising from the underlying grid, but actually sets it to the difference of the function value and the value in the given vector. This is commonly used in a time dependent, nonlinear problem for which we would like the rhs boundary values to be: U^{n+1}_k - U^{n+1}_{k+1} where n is the time iteration index, and k is the Newton iteration index. This means that the solution will be updated to U^{n+1}_{k+1} if the Jacobian is the identity for that row. This is very useful for time dependent boundary conditions for which the traditional method of letting the rhs value be zero does not work. Collective on GVec Input Parameters: + bd - The marker for the boundary to apply conditions along . u - A grid vector, usually the previous solution . field - The field to which the conditions apply . f - The function which defines the boundary condition . order - The test variable ordering - ctx - The user-supplied context Output Parameter: . b - The Rhs vector Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetVecBoundaryDifferenceRectangular(int bd, int field, GVec u, PointFunction f, VarOrdering order, GVec b, void *ctx) { Grid grid; Mesh mesh; int comp; /* The number of field components */ int size; /* The number of nodes in the boundary */ int *localStart; /* The offset of this field on a node of a given class */ int node; /* The canonical node number of the current boundary node */ int nclass; /* The class of the current boundary node */ double *x, *y, *z; /* Coordinates of the boundary nodes */ int vars; /* The number of variables affected (var/node * size) */ int *offsets; /* The canonical variable number for the first variable on each node */ int *rows; /* Rows corresponding to boundary nodes */ PetscScalar *values; /* Boundary values */ PetscScalar *uArray; /* The values in the vector u */ int firstVar; /* The canonical number of the first variable in this domain */ int rank; int i, j, count; #ifdef PETSC_USE_BOPT_g PetscTruth opt; #endif int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(b, VEC_COOKIE); PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE); ierr = GVecGetGrid(b, &grid); CHKERRQ(ierr); GridValidField(grid, field); ierr = MPI_Comm_rank(grid->comm, &rank); CHKERRQ(ierr); mesh = grid->mesh; comp = grid->fields[field].disc->comp; firstVar = order->firstVar[rank]; offsets = order->offsets; localStart = order->localStart[field]; /* Support for constrained problems */ ierr = VecGetSize(b, &size); CHKERRQ(ierr); if (grid->isConstrained) { if (size != grid->constraintOrder->numVars) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->constraintOrder->numVars); } offsets = grid->constraintOrder->offsets; } else { if (size != grid->order->numVars) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->order->numVars); } } /* Allocate memory */ ierr = GridGetBoundarySize(grid, bd, field, &size); CHKERRQ(ierr); if (size == 0) { #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscSynchronizedFlush(grid->comm); PetscSynchronizedFlush(grid->comm); } #endif ierr = VecAssemblyBegin(b); CHKERRQ(ierr); ierr = VecAssemblyEnd(b); CHKERRQ(ierr); PetscFunctionReturn(0); } vars = size*comp; ierr = PetscMalloc(size * sizeof(double), &x); CHKERRQ(ierr); ierr = PetscMalloc(size * sizeof(double), &y); CHKERRQ(ierr); ierr = PetscMalloc(size * sizeof(double), &z); CHKERRQ(ierr); ierr = PetscMalloc(vars * sizeof(PetscScalar), &values); CHKERRQ(ierr); ierr = PetscMalloc(vars * sizeof(int), &rows); CHKERRQ(ierr); /* Loop over boundary nodes */ ierr = GridGetBoundaryStart(grid, bd, field, PETSC_FALSE, &node, &nclass); CHKERRQ(ierr); for(i = 0, count = 0; node >= 0; i++) { for(j = 0; j < comp; j++, count++) { rows[count] = offsets[node] + j + localStart[nclass]; #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscSynchronizedPrintf(grid->comm, "[%d]bd %d field: %d node: %d row: %d class: %d\n", rank, bd, field, node, rows[count], nclass); } #endif } ierr = MeshGetNodeCoords(mesh, node, &x[i], &y[i], &z[i]); CHKERRQ(ierr); ierr = GridGetBoundaryNext(grid, bd, field, PETSC_FALSE, &node, &nclass); CHKERRQ(ierr); } #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscSynchronizedFlush(grid->comm); } #endif /* Get boundary values */ ierr = (*f)(size, comp, x, y, z, values, ctx); CHKERRQ(ierr); /* Taking the difference (we know that no values are off-processor) */ ierr = VecGetArray(u, &uArray); CHKERRQ(ierr); for(i = 0; i < vars; i++) values[i] = uArray[rows[i]-firstVar] - values[i]; ierr = VecRestoreArray(u, &uArray); CHKERRQ(ierr); /* Put values in Rhs */ #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscPrintf(grid->comm, "Setting boundary values in rhs bd %d field %d\n", bd, field); for(i = 0; i < vars; i++) PetscSynchronizedPrintf(grid->comm, " row: %d val: %g\n", rows[i], PetscRealPart(values[i])); PetscSynchronizedFlush(grid->comm); } ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif ierr = VecSetValues(b, vars, rows, values, INSERT_VALUES); CHKERRQ(ierr); ierr = VecAssemblyBegin(b); CHKERRQ(ierr); ierr = VecAssemblyEnd(b); CHKERRQ(ierr); ierr = PetscFree(x); CHKERRQ(ierr); ierr = PetscFree(y); CHKERRQ(ierr); ierr = PetscFree(z); CHKERRQ(ierr); ierr = PetscFree(values); CHKERRQ(ierr); ierr = PetscFree(rows); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetVecPointBoundaryDifference" /*@C GridSetVecPointBoundaryDifference This function sets Dirchlet boundary conditions on the linear Rhs arising from the underlying grid, but actually sets it to the difference of the function value and the value in the given vector. This is commonly used in a time dependent, nonlinear problem for which we would like the rhs boundary values to be: U^{n+1}_k - U^{n+1}_{k+1} where n is the time iteration index, and k is the Newton iteration index. This means that the solution will be updated to U^{n+1}_{k+1} if the Jacobian is the identity for that row. This is very useful for time dependent boundary conditions for which the traditional method of letting the rhs value be zero does not work. Collective on GVec Input Parameters: + node - The constrained node . u - A grid vector, usually the previous solution . field - The field to which the conditions apply . f - The function which defines the boundary condition - ctx - The user-supplied context Output Parameter: . b - The Rhs vector Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetVecPointBoundaryDifference(int node, int field, GVec u, PointFunction f, GVec b, void *ctx) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(b, VEC_COOKIE); ierr = GVecGetGrid(b, &grid); CHKERRQ(ierr); ierr = GridSetVecPointBoundaryDifferenceRectangular(node, field, u, f, grid->order, b, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridSetVecPointBoundaryDifferenceRectangular" /*@C GridSetVecBoundaryDifferenceRectangular This function sets Dirchlet boundary conditions on the linear Rhs arising from the underlying grid, but actually sets it to the difference of the function value and the value in the given vector. This is commonly used in a time dependent, nonlinear problem for which we would like the rhs boundary values to be: U^{n+1}_k - U^{n+1}_{k+1} where n is the time iteration index, and k is the Newton iteration index. This means that the solution will be updated to U^{n+1}_{k+1} if the Jacobian is the identity for that row. This is very useful for time dependent boundary conditions for which the traditional method of letting the rhs value be zero does not work. Collective on GVec Input Parameters: + node - The constrained node . u - A grid vector, usually the previous solution . field - The field to which the conditions apply . f - The function which defines the boundary condition . order - The test variable ordering - ctx - The user-supplied context Output Parameter: . b - The Rhs vector Level: advanced .keywords boundary conditions, finite element .seealso MeshGetBoundaryStart @*/ int GridSetVecPointBoundaryDifferenceRectangular(int node, int field, GVec u, PointFunction f, VarOrdering order, GVec b, void *ctx) { Grid grid; Mesh mesh; int comp; /* The number of field components */ int size; /* The number of nodes in the boundary */ int *localStart; /* The offset of this field on a node of a given class */ int nclass; /* The class of the current boundary node */ double x, y, z; /* Coordinates of the boundary nodes */ int *offsets; /* The canonical variable number for the first variable on each node */ int *rows; /* Rows corresponding to boundary nodes */ PetscScalar *values; /* Boundary values */ PetscScalar *uArray; /* The values in the vector u */ int firstVar; /* The canonical number of the first variable in this domain */ int rank; int i, j; #ifdef PETSC_USE_BOPT_g PetscTruth opt; #endif int ierr; PetscFunctionBegin; if (node < 0) { ierr = VecAssemblyBegin(b); CHKERRQ(ierr); ierr = VecAssemblyEnd(b); CHKERRQ(ierr); PetscFunctionReturn(0); } PetscValidHeaderSpecific(b, VEC_COOKIE); ierr = GVecGetGrid(b, &grid); CHKERRQ(ierr); ierr = GridGetMesh(grid, &mesh); CHKERRQ(ierr); GridValidField(grid, field); ierr = MPI_Comm_rank(grid->comm, &rank); CHKERRQ(ierr); comp = grid->fields[field].disc->comp; firstVar = order->firstVar[rank]; offsets = order->offsets; localStart = order->localStart[field]; /* Support for constrained problems */ ierr = VecGetSize(b, &size); CHKERRQ(ierr); if (grid->isConstrained) { if (size != grid->constraintOrder->numVars) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->constraintOrder->numVars); } offsets = grid->constraintOrder->offsets; } else { if (size != grid->order->numVars) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->order->numVars); } } /* Allocate memory */ size = 1; ierr = PetscMalloc(comp * sizeof(PetscScalar), &values); CHKERRQ(ierr); ierr = PetscMalloc(comp * sizeof(int), &rows); CHKERRQ(ierr); ierr = MeshGetNodeCoords(mesh, node, &x, &y, &z); CHKERRQ(ierr); ierr = GridGetNodeClass(grid, node, &nclass); CHKERRQ(ierr); for(j = 0; j < comp; j++) { rows[j] = offsets[node] + j + localStart[nclass]; #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscPrintf(PETSC_COMM_SELF, "[%d]field: %d node: %d row: %d class: %d\n", rank, field, node, rows[j], nclass); } #endif } /* Get boundary values */ ierr = (*f)(size, comp, &x, &y, &z, values, ctx); CHKERRQ(ierr); /* Taking the difference (we know that no values are off-processor) */ ierr = VecGetArray(u, &uArray); CHKERRQ(ierr); for(i = 0; i < comp; i++) values[i] = uArray[rows[i]-firstVar] - values[i]; ierr = VecRestoreArray(u, &uArray); CHKERRQ(ierr); /* Put values in Rhs */ #ifdef PETSC_USE_BOPT_g ierr = PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { PetscPrintf(grid->comm, "Setting boundary values on rhs node %d field %d\n", node, field); for(i = 0; i < comp; i++) PetscPrintf(PETSC_COMM_SELF, " row: %d val: %g\n", rows[i], PetscRealPart(values[i])); } ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__); CHKERRQ(ierr); #endif ierr = VecSetValues(b, comp, rows, values, INSERT_VALUES); CHKERRQ(ierr); ierr = VecAssemblyBegin(b); CHKERRQ(ierr); ierr = VecAssemblyEnd(b); CHKERRQ(ierr); ierr = PetscFree(values); CHKERRQ(ierr); ierr = PetscFree(rows); CHKERRQ(ierr); PetscFunctionReturn(0); }