Actual source code: gridBC.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: gridBC.c,v 1.3 2000/01/30 18:27:13 huangp Exp $";
  3: #endif

  5: #include "src/grid/gridimpl.h"     /*I "grid.h" I*//*I "gvec.h" I*/
  6: #include "src/vec/impls/mpi/pvecimpl.h" /* For GridCalcBCValues(), soon will not be needed */

  8: /*------------------------------------------------ Standard Functions -----------------------------------------------*/
  9: #undef  __FUNCT__
 11: /*@
 12:   GridDuplicateBC - Duplicates the boundary conditions of one grid in another.

 14:   Collective on Grid

 16:   Input Parameter:
 17: . grid    - The grid

 19:   Output Parameter:
 20: . newGrid - The altered grid

 22:   Level: intermediate

 24: .keywords: grid, duplicate, BC
 25: .seealso: GridDuplicate()
 26: @*/
 27: int GridDuplicateBC(Grid grid, Grid newGrid)
 28: {
 29:   int bc;

 35:   for(bc = 0; bc < grid->numBC; bc++) {
 36:     GridAddBC(newGrid, grid->bc[bc].boundary, grid->bc[bc].field, grid->bc[bc].func, grid->bc[bc].reduce);
 37:   }
 38:   for(bc = 0; bc < grid->numPointBC; bc++) {
 39:     GridAddPointBC(newGrid, grid->bc[bc].point[0], grid->bc[bc].point[1], grid->bc[bc].point[2], grid->bc[bc].field,
 40:                           grid->bc[bc].func, grid->bc[bc].reduce);
 41:   }
 42:   return(0);
 43: }

 45: #undef  __FUNCT__
 47: /*@
 48:   GridFinalizeBC - Destroys all structures associated with explicit system
 49:   reduction using boundary conditions. This should be called after all
 50:   calculation is finished, prior to GridDestroy().

 52:   Collective on Grid

 54:   Input Parameter:
 55: . grid - The grid

 57:   Level: beginner

 59: .keywords: grid, boundary conditions
 60: .seealso: GridSetBC(), GridAddBC()
 61: @*/
 62: int GridFinalizeBC(Grid grid)
 63: {

 68:   if (grid->bdReduceVec) {
 69:     GVecDestroy(grid->bdReduceVec);
 70:   }
 71:   if (grid->bdReduceVecOld) {
 72:     GVecDestroy(grid->bdReduceVecOld);
 73:   }
 74:   if (grid->bdReduceVecDiff) {
 75:     GVecDestroy(grid->bdReduceVecDiff);
 76:   }
 77:   if (grid->bdReduceMat) {
 78:     GMatDestroy(grid->bdReduceMat);
 79:   }
 80:   if (grid->reduceVec) {
 81:     GVecDestroy(grid->reduceVec);
 82:   }

 84:   return(0);
 85: }

 87: /*----------------------------------------------- Database Functions ------------------------------------------------*/
 88: #undef  __FUNCT__
 90: /*@C GridSetBC
 91:   This function sets the boundary condition to use for the problem.

 93:   Collective on Grid

 95:   Input Parameters:
 96: + grid   - The grid
 97: . bd     - The marker for the boundary along which conditions are applied
 98: . field  - The field to which the boundary condition is applied
 99: . f      - The function which defines the boundary condition
100: - reduce - The flag for explicit reduction of the system

102:   Level: intermediate

104: .keywords active field
105: .seealso GridAddBC, GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction
106: @*/
107: int GridSetBC(Grid grid, int bd, int field, PointFunction f, PetscTruth reduce)
108: {

113:   GridValidField(grid, field);
114:   grid->numBC = 0;
115:   GridAddBC(grid, bd, field, f, reduce);
116:   return(0);
117: }

119: #undef  __FUNCT__
121: /*@C GridAddBC
122:   This function adds a boundary condition to use for the problem.

124:   Collective on Grid

126:   Input Parameters:
127: + grid   - The grid
128: . bd     - The marker for the boundary along which conditions are applied
129: . field  - The field to which the boundary condition is applied
130: . f      - The function which defines the boundary condition
131: - reduce - The flag for explicit reduction of the system

133:   Level: intermediate

135: .keywords active field
136: .seealso GridSetBC, GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction
137: @*/
138: int GridAddBC(Grid grid, int bd, int field, PointFunction f, PetscTruth reduce)
139: {
140:   GridBC *tempBC;
141:   int     bdIndex;
142:   int     ierr;

146:   GridValidField(grid, field);
147:   while (grid->numBC >= grid->maxBC) {
148:     PetscMalloc(grid->maxBC*2 * sizeof(GridBC), &tempBC);
149:     PetscLogObjectMemory(grid, grid->maxBC * sizeof(GridBC));
150:     PetscMemcpy(tempBC, grid->bc, grid->maxBC * sizeof(GridBC));
151:     PetscFree(grid->bc);
152:     grid->bc     = tempBC;
153:     grid->maxBC *= 2;
154:   }
155:   /* Make sure boundary is legal */
156:   MeshGetBoundaryIndex(grid->mesh, bd, &bdIndex);
157:   grid->bc[grid->numBC].boundary = bd;
158:   grid->bc[grid->numBC].field    = field;
159:   grid->bc[grid->numBC].func     = f;
160:   grid->bc[grid->numBC].reduce   = reduce;
161:   grid->bc[grid->numBC].node     = -1;
162:   grid->numBC++;
163:   /* Check whether to reduce system */
164:   if (reduce == PETSC_TRUE) grid->reduceSystem = PETSC_TRUE;
165:   return(0);
166: }

168: #undef  __FUNCT__
170: /*@C GridSetPointBC
171:   This function sets the boundary condition to use for the problem at a point.

173:   Collective on Grid

175:   Input Parameters:
176: + grid   - The grid
177: . x,y,z  - The point at which conditions are applied
178: . field  - The field to which the boundary condition is applied
179: . f      - The function which defines the boundary condition
180: - reduce - The flag for explicit reduction of the system

182:   Level: intermediate

184: .keywords active field
185: .seealso GridAddBC, GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction
186: @*/
187: int GridSetPointBC(Grid grid, double x, double y, double z, int field, PointFunction f, PetscTruth reduce)
188: {

193:   grid->numPointBC = 0;
194:   GridAddPointBC(grid, x, y, z, field, f, reduce);
195:   return(0);
196: }

198: #undef  __FUNCT__
200: /*@C GridAddPointBC
201:   This function adds a boundary condition to use for the problem at a point.

203:   Collective on Grid

205:   Input Parameters:
206: + grid   - The grid
207: . x,y,z  - The point at which conditions are applied
208: . field  - The field to which the boundary condition is applied
209: . f      - The function which defines the boundary condition
210: - reduce - The flag for explicit reduction of the system

212:   Level: intermediate

214: .keywords active field
215: .seealso GridSetBC, GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction
216: @*/
217: int GridAddPointBC(Grid grid, double x, double y, double z, int field, PointFunction f, PetscTruth reduce)
218: {
219:   GridBC *tempBC;
220:   int     ierr;

224:   GridValidField(grid, field);
225:   while (grid->numPointBC >= grid->maxPointBC) {
226:     PetscMalloc(grid->maxPointBC*2 * sizeof(GridBC), &tempBC);
227:     PetscLogObjectMemory(grid, grid->maxPointBC * sizeof(GridBC));
228:     PetscMemcpy(tempBC, grid->pointBC, grid->maxPointBC * sizeof(GridBC));
229:     PetscFree(grid->pointBC);
230:     grid->pointBC     = tempBC;
231:     grid->maxPointBC *= 2;
232:   }
233:   if (GridGetNearestBdNode(grid, field, x, y, z, &grid->pointBC[grid->numPointBC].node)) {
234:     SETERRQ3(PETSC_ERR_ARG_WRONG, "Invalid point {%g,%g,%g} specified for boundary condition", x, y, z);
235:   }
236:   grid->pointBC[grid->numPointBC].point[0] = x;
237:   grid->pointBC[grid->numPointBC].point[1] = y;
238:   grid->pointBC[grid->numPointBC].point[2] = z;
239:   grid->pointBC[grid->numPointBC].field    = field;
240:   grid->pointBC[grid->numPointBC].func     = f;
241:   grid->pointBC[grid->numPointBC].reduce   = reduce;
242:   grid->pointBC[grid->numPointBC].boundary = -1;
243:   grid->numPointBC++;
244:   /* Check whether to reduce system */
245:   if (reduce == PETSC_TRUE) grid->reduceSystem = PETSC_TRUE;
246:   return(0);
247: }

249: #undef  __FUNCT__
251: /*@
252:   GridSetBCMultiplier - This sets the scalar multiplier used for reduction components on the rhs.

254:   Collective on Grid

256:   Input Parameters:
257: + grid  - The grid
258: - alpha - The scalar multiplier

260:   Note:
261:   For example, this should be -1 in a nonlinear iteration. The default is 1.

263:   Level: developer

265: .keywords: grid, reduction, boundary conditions
266: .seealso: GridGetBCMultiplier(), GridSetBC(), GridAddBC()
267: @*/
268: int GridSetBCMultiplier(Grid grid, PetscScalar alpha)
269: {
272:   grid->reduceAlpha = alpha;
273:   return(0);
274: }

276: #undef  __FUNCT__
278: /*@
279:   GridGetBCMultiplier - This gets the scalar multiplier used for reduction components on the rhs.

281:   Not collective

283:   Input Parameter:
284: . grid  - The grid

286:   Output Parameter:
287: . alpha - The scalar multiplier

289:   Level: developer

291: .keywords: grid, reduction, boundary conditions
292: .seealso: GridSetBCMultiplier(), GridSetBC(), GridAddBC()
293: @*/
294: int GridGetBCMultiplier(Grid grid, PetscScalar *alpha)
295: {
299:   *alpha = grid->reduceAlpha;
300:   return(0);
301: }

303: #undef  __FUNCT__
305: /*@
306:   GridSetBCContext - This sets the optional user context passed to all
307:   routines which assemble boundary reduction information. Must be called
308:   before GridSetUp().

310:   Collective on Grid

312:   Input Parameters:
313: + grid - The grid
314: - ctx  - The context

316:   Level: intermediate

318: .keywords: grid, reduction, boundary conditions
319: .seealso: GridGetBCContext(), GridSetBC(), GridAddBC()
320: @*/
321: int GridSetBCContext(Grid grid, void *ctx)
322: {
325:   grid->reduceContext = ctx;
326:   return(0);
327: }

329: #undef  __FUNCT__
331: /*@
332:   GridGetBCContext - This gets the optional user context passed to all
333:   routines which assemble boundary reduction information.

335:   Not collective

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

340:   Output parameter:
341: . ctx  - The context

343:   Level: intermediate

345: .keywords: grid, reduction, boundary conditions
346: .seealso: GridSetBCContext(), GridSetBC(), GridAddBC()
347: @*/
348: int GridGetBCContext(Grid grid, void **ctx)
349: {
353:   *ctx = grid->reduceContext;
354:   return(0);
355: }

357: #undef  __FUNCT__
359: /*@
360:   GridSetBCValuesType - This determines which boundary values are used to reduce
361:   the system. It is intended to allow time dependent boundary conditions to be
362:   used, and also supports the difference of two sets of values.

364:   Collective on Grid

366:   Input Parameter:
367: . grid - The grid

369:   Level: intermediate

371: .keywords: grid, reduction, boundary conditions
372: .seealso: GridSetBC(), GridAddBC()
373: @*/
374: int GridSetBCValuesType(Grid grid, BCValuesType type)
375: {
378:   if (grid->reduceSystem == PETSC_FALSE)
379:     return(0);

381:   switch(type) {
382:   case BC_VALUES:
383:     grid->bdReduceVecCur = grid->bdReduceVec;
384:     break;
385:   case BC_VALUES_OLD:
386:     if (grid->bdReduceVecOld == PETSC_NULL) {
387:       SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Old boundary values not stored");
388:     }
389:     grid->bdReduceVecCur = grid->bdReduceVecOld;
390:     break;
391:   case BC_VALUES_DIFF:
392:     if (grid->bdReduceVecDiff == PETSC_NULL) {
393:       SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Difference of boundary values not stored");
394:     }
395:     grid->bdReduceVecCur = grid->bdReduceVecDiff;
396:     break;
397:   default:
398:     SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid type %d for boundary value calculation", type);
399:   }
400:   return(0);
401: }

403: /*---------------------------------------------- Calculation Functions ----------------------------------------------*/
404: #undef  __FUNCT__
406: /*@C GridCalcPointBCNodes
407:   This function recalculates the nodes used for point boundary conditions.

409:   Collective on Grid

411:   Input Parameter:
412: . grid   - The grid

414:   Notes:
415:   This function is called by GridReform() after the mesh is recalculated.

417:   Level: advanced

419: .keywords grid, point BC, node
420: .seealso GridSetBC, GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction
421: @*/
422: int GridCalcPointBCNodes(Grid grid)
423: {
424:   double x, y, z;
425:   int    bc;

428:   for(bc = 0; bc < grid->numPointBC; bc++) {
429:     x = grid->pointBC[bc].point[0];
430:     y = grid->pointBC[bc].point[1];
431:     z = grid->pointBC[bc].point[2];
432:     if (GridGetNearestBdNode(grid, grid->pointBC[bc].field, x, y, z, &grid->pointBC[bc].node)) {
433:       SETERRQ3(PETSC_ERR_ARG_WRONG, "Invalid point {%g,%g,%g} specified for boundary condition", x, y, z);
434:     }
435:   }
436:   return(0);
437: }

439: #undef  __FUNCT__
441: int GridSaveBCValues_Private(Grid grid, VarOrdering reduceOrder, Vec reduceVec) {
442:   PetscScalar *array, *arrayOld;
443:   int          ierr;

446:   /* Create storage for reduction of Rhs */
447:   if (grid->bdReduceVecOld == PETSC_NULL) {
448:     GVecDuplicate(reduceVec, &grid->bdReduceVecOld);
449:   } else if (((Vec_MPI *) grid->bdReduceVecOld->data)->nghost != ((Vec_MPI *) grid->bdReduceVec->data)->nghost) {
450:     GVecDestroy(grid->bdReduceVecOld);
451:     GVecDuplicate(reduceVec, &grid->bdReduceVecOld);
452:   }
453:   /* VecCopy(grid->bdReduceVec, grid->bdReduceVecOld);                                              */
454:   VecGetArray(reduceVec,            &array);
455:   VecGetArray(grid->bdReduceVecOld, &arrayOld);
456:   PetscMemcpy(arrayOld, array, reduceOrder->numOverlapVars * sizeof(PetscScalar));
457:   VecRestoreArray(reduceVec,            &array);
458:   VecRestoreArray(grid->bdReduceVecOld, &arrayOld);
459:   return(0);
460: }

462: #undef  __FUNCT__
464: int GridCalcGridBCValues_Private(Grid grid, VarOrdering reduceOrder, Vec reduceVec, void *ctx) {
465:   GridBC     *gBC = grid->bc;
466:   VarOrdering bcOrder;
467:   int         bc;
468:   int         ierr;

471:   /* Evaluate the vector of boundary values --
472:        If order->localStart[field] is NULL, this means the field is not present in the ordering. This is
473:      a better check than seeing if the field is active, since we might want to pass in an order on that
474:      field to make boundary values for an inactive field.
475:   */
476:   for(bc = 0; bc < grid->numBC; bc++) {
477:     if (gBC[bc].reduce                         != PETSC_TRUE) continue;
478:     if (reduceOrder->localStart[gBC[bc].field] == PETSC_NULL) continue;
479:     VarOrderingCreateSubset(reduceOrder, 1, &gBC[bc].field, PETSC_FALSE, &bcOrder);
480:     (*grid->ops->gvecevaluatefunctionboundary)(grid, reduceVec, gBC[bc].boundary, bcOrder, gBC[bc].func, 1.0, ctx);
481:     VarOrderingDestroy(bcOrder);
482: #ifdef PETSC_USE_BOPT_g
483:     PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
484: #endif
485:   }
486:   return(0);
487: }

489: #undef  __FUNCT__
491: int GridCalcPointBCValues_Private(Grid grid, VarOrdering reduceOrder, Vec reduceVec, void *ctx) {
492:   GridBC       *pBC          = grid->pointBC;
493:   int         **localStart   = reduceOrder->localStart;
494:   int          *offsets      = reduceOrder->offsets;
495:   int          *localOffsets = reduceOrder->localOffsets;
496:   FieldClassMap map;
497:   VarOrdering   bcOrder;
498:   PetscScalar  *array;
499:   double        x, y, z;
500:   int           numNodes, firstVar, rank;
501:   int           bc, field, numComp, node, nclass, row;
502:   int           ierr;

505:   MPI_Comm_rank(grid->comm, &rank);
506:   VarOrderingGetClassMap(reduceOrder, &map);
507:   numNodes = map->numNodes;
508:   firstVar = reduceOrder->firstVar[rank];
509:   /* Evaluate the vector of boundary values --
510:        If order->localStart[field] is NULL, this means the field is not present in the ordering. This is
511:      a better check than seeing if the field is active, since we might want to pass in an order on that
512:      field to make boundary values for an inactive field.
513:   */
514:   VecGetArray(reduceVec, &array);
515:   for(bc = 0; bc < grid->numPointBC; bc++) {
516:     /* BC is not reduced out of the system */
517:     if (pBC[bc].reduce                         != PETSC_TRUE) continue;
518:     /* Field is not present in the ordering */
519:     if (reduceOrder->localStart[pBC[bc].field] == PETSC_NULL) continue;

521:     field   = pBC[bc].field;
522:     numComp = grid->fields[field].numComp;
523:     VarOrderingCreateSubset(reduceOrder, 1, &field, PETSC_FALSE, &bcOrder);

525:     /* Point is in another domain */
526:     if (pBC[bc].node < 0) continue;
527:     node   = pBC[bc].node;
528:     nclass = map->classes[node];

530:     if (node >= numNodes) {
531:       row = localOffsets[node-numNodes];
532:     } else {
533:       row = offsets[node] - firstVar + localStart[field][nclass];
534:     }
535:     x = pBC[bc].point[0];
536:     y = pBC[bc].point[1];
537:     z = pBC[bc].point[2];
538:     (*pBC[bc].func)(1, numComp, &x, &y, &z, &array[row], ctx);

540:     VarOrderingDestroy(bcOrder);
541: #ifdef PETSC_USE_BOPT_g
542:     PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
543: #endif
544:   }
545:   VecRestoreArray(reduceVec, &array);
546:   return(0);
547: }

549: #undef  __FUNCT__
551: int GridCalcBCValues_Private(Grid grid, VarOrdering reduceOrder, Vec reduceVec, PetscTruth save, void *ctx) {
552:   PetscScalar *array;
553:   int          numGhostVars;
554:   int          ierr;

557:   numGhostVars = reduceOrder->numOverlapVars - reduceOrder->numLocVars;
558:   if (((Vec_MPI *) reduceVec->data)->nghost != numGhostVars) {
559:     SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid reduce vector size %d should be %d", ((Vec_MPI *) reduceVec->data)->nghost, numGhostVars);
560:   }

562:   if (save == PETSC_TRUE) {
563:     GridSaveBCValues_Private(grid, reduceOrder, reduceVec);
564:   }

566:   /* Initialize vector */
567:   /* VecSet(&zero, reduceVec);                                                                      */
568:   VecGetArray(reduceVec, &array);
569:   PetscMemzero(array, reduceOrder->numOverlapVars * sizeof(PetscScalar));
570:   VecRestoreArray(reduceVec, &array);

572:   GridCalcGridBCValues_Private(grid, reduceOrder, reduceVec, ctx);
573:   GridCalcPointBCValues_Private(grid, reduceOrder, reduceVec, ctx);

575:   return(0);
576: }

578: #undef  __FUNCT__
580: /*@
581:   GridCalcBCValues - This function calculates the boundary values. It
582:   is normally called once a timestep when using time dependent boundary
583:   conditions.

585:   Collective on Grid

587:   Input Parameters:
588: + grid - The grid
589: . save - A flag used to store old values, usually for timestepping
590: - ctx  - The context

592:   Level: advanced

594: .keywords: grid, reduction, boundary conditions
595: .seealso: GridSetBCContext(), GridSetBC(), GridAddBC()
596: @*/
597: int GridCalcBCValues(Grid grid, PetscTruth save, void *ctx)
598: {


604:   if (grid->reduceSystem == PETSC_TRUE) {
606:     GridCalcBCValues_Private(grid, grid->reduceOrder, grid->bdReduceVec, save, ctx);
607:   }

609:   return(0);
610: }

612: #undef  __FUNCT__
614: /*@
615:   GridCalcBCValuesDifference - This function calculates the difference of the
616:   last two sets of boundary values and puts it in an internal vector. This is
617:   is normally used to implement the Rhs time derivative in a GTS.

619:   Collective on Grid

621:   Input Parameter:
622: . grid - The grid

624:   Level: advanced

626: .keywords: grid, reduction, boundary conditions
627: .seealso: GridSetBCContext(), GridSetBC(), GridAddBC()
628: @*/
629: int GridCalcBCValuesDifference(Grid grid)
630: {
631:   PetscScalar *array, *arrayOld, *arrayDiff;
632:   int          numGhostVars;
633:   register int i, n;
634:   int          ierr;

638:   if (grid->reduceSystem == PETSC_TRUE) {
639:     numGhostVars = grid->reduceOrder->numOverlapVars - grid->reduceOrder->numLocVars;
640:     if (((Vec_MPI *) grid->bdReduceVec->data)->nghost != numGhostVars) {
641:       SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduce vector size");
642:     }
643:     if (grid->bdReduceVecOld == PETSC_NULL) {
644:       SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "No previous boundary values");
645:     }
646:     /* Create storage for reduction of Rhs */
647:     if (grid->bdReduceVecDiff == PETSC_NULL) {
648:       GVecDuplicate(grid->bdReduceVec, &grid->bdReduceVecDiff);
649:     } else if (((Vec_MPI *) grid->bdReduceVecDiff->data)->nghost != ((Vec_MPI *) grid->bdReduceVec->data)->nghost) {
650:       GVecDestroy(grid->bdReduceVecDiff);
651:       GVecDuplicate(grid->bdReduceVec, &grid->bdReduceVecDiff);
652:     }
653:     /* VecWAXPY(&minusOne, grid->bdReduceVecOld, grid->bdReduceVec, grid->bdReduceVecDiff);           */
654:     VecGetArray(grid->bdReduceVec,     &array);
655:     VecGetArray(grid->bdReduceVecOld,  &arrayOld);
656:     VecGetArray(grid->bdReduceVecDiff, &arrayDiff);
657:     n    = grid->reduceOrder->numOverlapVars;
658:     PetscLogFlops(n);
659:     for(i = 0; i < n; i++)
660:       arrayDiff[i] = array[i] - arrayOld[i];
661:     VecRestoreArray(grid->bdReduceVec,     &array);
662:     VecRestoreArray(grid->bdReduceVecOld,  &arrayOld);
663:     VecRestoreArray(grid->bdReduceVecDiff, &arrayDiff);
664:   }
665: #ifdef PETSC_USE_BOPT_g
666:   PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
667: #endif

669:   return(0);
670: }

672: /*----------------------------------------------- Reduction Functions -----------------------------------------------*/
673: #undef  __FUNCT__
675: /*@C
676:   GridSetReduceSystem - This function determines whether unknowns associated
677:   with boundary conditions are eliminated from the system.

679:   Collective on Grid

681:   Input Parameters:
682: + grid   - The grid
683: - reduce - The flag for explicit reduction of the system

685:   Level: intermediate

687: .keywords grid, boundary condition, reduce
688: .seealso GridGetReduceSystem(), GridSetBC(), GridAddBC(), GridSetPointBC(), GridAddPointBC()
689: @*/
690: int GridSetReduceSystem(Grid grid, PetscTruth reduce)
691: {
694:   grid->reduceSystem = reduce;
695:   return(0);
696: }

698: #undef  __FUNCT__
700: /*@C
701:   GridGetReduceSystem - This function reveals whether unknowns associated
702:   with boundary conditions are eliminated from the system.

704:   Not collective

706:   Input Parameter:
707: . grid   - The grid

709:   Output Parameter:
710: . reduce - The flag for explicit reduction of the system

712:   Level: intermediate

714: .keywords grid, boundary condition, reduce
715: .seealso GridSetReduceSystem(), GridSetBC(), GridAddBC(), GridSetPointBC(), GridAddPointBC()
716: @*/
717: int GridGetReduceSystem(Grid grid, PetscTruth *reduce)
718: {
722:   *reduce = grid->reduceSystem;
723:   return(0);
724: }

726: #undef  __FUNCT__
728: /*@C
729:   GridSetReduceElement - This function determines whether element matrices and vectors
730:   are reduced on the fly, or if boundary operators are stored and applied.

732:   Collective on Grid

734:   Input Parameters:
735: + grid   - The grid
736: - reduce - The flag for explicit reduction of the system

738:   Level: intermediate

740: .keywords grid, boundary condition, reduce, element
741: .seealso GridGetReduceElement(), GridSetBC(), GridAddBC(), GridSetPointBC(), GridAddPointBC()
742: @*/
743: int GridSetReduceElement(Grid grid, PetscTruth reduce)
744: {
747:   grid->reduceElement = reduce;
748:   return(0);
749: }

751: #undef  __FUNCT__
753: /*@C
754:   GridGetReduceElement - This function indicates whether element matrices and vectors
755:   are reduced on the fly, or if boundary operators are stored and applied.

757:   Not collective

759:   Input Parameter:
760: . grid   - The grid

762:   Output Parameter:
763: . reduce - The flag for explicit reduction of the system

765:   Level: intermediate

767: .keywords grid, boundary condition, reduce, element
768: .seealso GridSetReduceElement(), GridSetBC(), GridAddBC(), GridSetPointBC(), GridAddPointBC()
769: @*/
770: int GridGetReduceElement(Grid grid, PetscTruth *reduce)
771: {
775:   *reduce = grid->reduceElement;
776:   return(0);
777: }

779: /*---------------------------------------------- Application Functions ----------------------------------------------*/
780: #undef  __FUNCT__
782: /*
783:   GridResetConstrainedMultiply_Private - This function resets the mulplication routine for constrained matrices

785:   Input Parameters:
786: + grid - The Grid
787: - A    - The GMat

789:   Level: developer

791: .keywords Grid, GMat, reset, constrained, multiply
792: .seealso GridEvaluateRhs
793: */
794: int GridResetConstrainedMultiply_Private(Grid grid, GMat A) {
795:   PetscTruth isConstrained, explicitConstraints;
796:   void     (*oldMult)(void);
797:   int        ierr;

800:   GridIsConstrained(grid, &isConstrained);
801:   GridGetExplicitConstraints(grid, &explicitConstraints);
802:   if (isConstrained == PETSC_TRUE) {
803:     if (explicitConstraints == PETSC_FALSE) {
804:       MatShellGetOperation(A, MATOP_MULT, &oldMult);
805:       if (oldMult != (void (*)(void)) GMatMatMultConstrained) {
806:         MatShellSetOperation(A, MATOP_MULT_CONSTRAINED, oldMult);
807:       }
808:       MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) GMatMatMultConstrained);

810:       MatShellGetOperation(A, MATOP_MULT_TRANSPOSE, &oldMult);
811:       if (oldMult != (void (*)(void)) GMatMatMultTransposeConstrained) {
812:         MatShellSetOperation(A, MATOP_MULT_TRANSPOSE_CONSTRAINED, oldMult);
813:       }
814:       MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, (void (*)(void)) GMatMatMultTransposeConstrained);
815:     } else {
816:       MatShellGetOperation(A, MATOP_MULT_CONSTRAINED, &oldMult);
817:       if (oldMult != PETSC_NULL) {
818:         MatShellSetOperation(A, MATOP_MULT, oldMult);
819:       }

821:       MatShellGetOperation(A, MATOP_MULT_TRANSPOSE_CONSTRAINED, &oldMult);
822:       if (oldMult != PETSC_NULL) {
823:         MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, oldMult);
824:       }
825:     }
826:   }
827:   return(0);
828: }

830: #undef  __FUNCT__
832: /*@C GridSetBoundary
833:   This function sets Dirchlet boundary conditions on the linear problem arising
834:   from the underlying grid.

836:   Collective on Grid

838:   Input Parameters:
839: + bd       - The marker for the boundary to apply conditions along
840: . field    - The field to which the conditions apply
841: . diag     - The scalar to be placed on the diagonal
842: . f        - The function which defines the boundary condition
843: - ctx      - The user-supplied context

845:   Output Parameters:
846: + A        - The system matrix
847: - b        - The Rhs vector

849:   Level: intermediate

851: .keywords boundary conditions, finite element
852: .seealso MeshGetBoundaryStart
853: @*/
854: int GridSetBoundary(int bd, int field, PetscScalar diag, PointFunction f, GMat A, GVec b, void *ctx)
855: {
856:   Grid grid;
857:   int  ierr;

861:   GMatGetGrid(A, &grid);
862:   GridSetBoundaryRectangular(bd, field, diag, f, grid->order, A, b, ctx);
863:   return(0);
864: }

866: #undef  __FUNCT__
868: /*@C GridSetBoundaryRectangular
869:   This function sets Dirchlet boundary conditions on the linear problem arising
870:   from the underlying grid, and the default variable ordering can be overridden.

872:   Collective on Grid

874:   Input Parameters:
875: + bd       - The marker for the boundary to apply conditions along
876: . field    - The field to which the conditions apply
877: . diag     - The scalar to be placed on the diagonal
878: . f        - The function which defines the boundary condition
879: . order    - The test variable ordering
880: - ctx      - The user-supplied context

882:   Output Parameters:
883: + A        - The system matrix
884: - b        - The Rhs vector

886:   Level: advanced

888: .keywords boundary conditions, finite element
889: .seealso MeshGetBoundaryStart
890: @*/
891: int GridSetBoundaryRectangular(int bd, int field, PetscScalar diag, PointFunction f, VarOrdering order, GMat A, GVec b, void *ctx)
892: {
893:   Grid           grid, grid2;
894:   Mesh           mesh;
895:   int            comp;       /* The number of field components */
896:   int            size;       /* The number of nodes in the boundary */
897:   int           *localStart; /* The offset of this field on a node of a given class */
898:   int            node;       /* The canonical node number of the current boundary node */
899:   int            nclass;     /* The class of the current boundary node */
900:   double        *x, *y, *z;  /* Coordinates of the boundary nodes */
901:   int            vars;       /* The number of variables affected (var/node * size) */
902:   int           *offsets;    /* The canonical variable number for the first variable on each node */
903:   int           *rows;       /* Rows corresponding to boundary nodes */
904:   PetscScalar   *values;     /* Boundary values */
905:   PetscScalar    elem = diag;
906:   IS             is;
907:   int            rank;
908:   int            i, j, count;
909: #ifdef PETSC_USE_BOPT_g
910:   PetscTruth     opt;
911: #endif
912:   int            ierr;

918:   GMatGetGrid(A, &grid);
919:   GVecGetGrid(b, &grid2);
920:   if (grid != grid2) SETERRQ(PETSC_ERR_ARG_INCOMP, "Matrix and vector have different underlying grids");
921:   GridValidField(grid, field);
922:   MPI_Comm_rank(grid->comm, &rank);
923:   mesh       = grid->mesh;
924:   comp       = grid->fields[field].disc->comp;
925:   offsets    = order->offsets;
926:   localStart = order->localStart[field];

928:   /* Support for constrained problems */
929:   VecGetSize(b, &size);
930:   if (grid->isConstrained) {
931:     if (size != grid->constraintOrder->numVars) SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid vector size");
932:     offsets = grid->constraintOrder->offsets;
933:   } else {
934:     if (size != grid->order->numVars) SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid vector size");
935:   }

937:   /* Allocate memory */
938:   GridGetBoundarySize(grid, bd, field, &size);
939:   if (size == 0) {
940: #ifdef PETSC_USE_BOPT_g
941:     PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
942:     if (opt == PETSC_TRUE) {
943:       PetscSynchronizedFlush(grid->comm);
944:       PetscSynchronizedFlush(grid->comm);
945:     }
946: #endif
947:     VecAssemblyBegin(b);
948:     VecAssemblyEnd(b);
949:     ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is);
950:     MatZeroRows(A, is, &elem);
951:     ISDestroy(is);
952:     return(0);
953:   }
954:   vars = size*comp;
955:   PetscMalloc(size * sizeof(double),      &x);
956:   PetscMalloc(size * sizeof(double),      &y);
957:   PetscMalloc(size * sizeof(double),      &z);
958:   PetscMalloc(vars * sizeof(PetscScalar), &values);
959:   PetscMalloc(vars * sizeof(int),         &rows);

961:   /* Loop over boundary nodes */
962:   GridGetBoundaryStart(grid, bd, field, PETSC_FALSE, &node, &nclass);
963:   for(i = 0, count = 0; node >= 0; i++) {
964:     for(j = 0; j < comp; j++, count++) {
965:       rows[count] = offsets[node] + j + localStart[nclass];
966: #ifdef PETSC_USE_BOPT_g
967:     PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
968:     if (opt == PETSC_TRUE) {
969:       PetscSynchronizedPrintf(grid->comm, "[%d]bd %d field: %d node: %d row: %d class: %dn",
970:                               rank, bd, field, node, rows[count], nclass);
971:     }
972: #endif
973:     }
974:     MeshGetNodeCoords(mesh, node, &x[i], &y[i], &z[i]);
975:     GridGetBoundaryNext(grid, bd, field, PETSC_FALSE, &node, &nclass);
976:   }
977: #ifdef PETSC_USE_BOPT_g
978:   PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
979:   if (opt == PETSC_TRUE) {
980:     PetscSynchronizedFlush(grid->comm);
981:   }
982: #endif
983:   /* Get boundary values */
984:   (*f)(size, comp, x, y, z, values, ctx);
985:   /* Put values in Rhs */
986: #ifdef PETSC_USE_BOPT_g
987:   PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
988:   if (opt == PETSC_TRUE) {
989:     PetscPrintf(grid->comm, "Setting boundary values in rhs bd %d field %dn", bd, field);
990:     for(i = 0; i < vars; i++) PetscSynchronizedPrintf(grid->comm, "  row: %d val: %gn", rows[i], PetscRealPart(values[i]));
991:     PetscSynchronizedFlush(grid->comm);
992:   }
993:   PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
994: #endif
995:   VecSetValues(b, vars, rows, values, INSERT_VALUES);
996:   VecAssemblyBegin(b);
997:   VecAssemblyEnd(b);
998:   /* Set rows of A to the identity */
999:   ISCreateGeneral(PETSC_COMM_SELF, vars, rows, &is);
1000:   MatZeroRows(A, is, &elem);
1001:   ISDestroy(is);

1003:   GridResetConstrainedMultiply_Private(grid, A);

1005:   PetscFree(x);
1006:   PetscFree(y);
1007:   PetscFree(z);
1008:   PetscFree(values);
1009:   PetscFree(rows);
1010:   return(0);
1011: }

1013: /*------------------------------------------------- Matrix Functions ------------------------------------------------*/
1014: #undef  __FUNCT__
1016: /*@C GridSetMatBoundary
1017:   This function sets Dirchlet boundary conditions on the linear system matrix arising
1018:   from the underlying grid.

1020:   Collective on GMat

1022:   Input Parameters:
1023: + bd       - The marker for the boundary to apply conditions along
1024: . field    - The field to which the conditions apply
1025: . diag     - The scalar to be placed on the diagonal
1026: - ctx      - The user-supplied context

1028:   Output Parameter:
1029: . A     - The system matrix

1031:   Level: advanced

1033: .keywords boundary conditions, finite element
1034: .seealso MeshGetBoundaryStart
1035: @*/
1036: int GridSetMatBoundary(int bd, int field, PetscScalar diag, GMat A, void *ctx)
1037: {
1038:   Grid grid;
1039:   int  ierr;

1043:   GMatGetGrid(A, &grid);
1044:   GridSetMatBoundaryRectangular(1, &bd, &field, diag, grid->order, A, ctx);
1045:   return(0);
1046: }

1048: #undef  __FUNCT__
1050: /*@C GridSetMatBoundaryRectangular
1051:   This function sets Dirchlet boundary conditions on the linear system matrix arising
1052:   from the underlying grid, and the default variable ordering can be overridden.

1054:   Collective on GMat

1056:   Input Parameters:
1057: + num   - The number of boundary conditions
1058: . bd    - The markers for each boundary to apply conditions along
1059: . field - The fields to which the conditions apply
1060: . diag  - The scalar to be placed on the diagonal
1061: . order - The test variable ordering
1062: - ctx   - The user-supplied context

1064:   Output Parameter:
1065: . A     - The system matrix

1067:   Level: advanced

1069: .keywords boundary conditions, finite element
1070: .seealso MeshGetBoundaryStart
1071: @*/
1072: int GridSetMatBoundaryRectangular(int num, int *bd, int *field, PetscScalar diag, VarOrdering order, GMat A, void *ctx)
1073: {
1074:   Grid           grid;
1075:   int            comp;       /* The number of field components */
1076:   int            size;       /* The number of nodes in the boundary */
1077:   int            totSize;    /* The number of nodes in all boundaries */
1078:   int           *localStart; /* The offset of this field on a node of a given class */
1079:   int            node;       /* The canonical node number of the current boundary node */
1080:   int            nclass;     /* The class of the current boundary node */
1081:   int            vars;       /* The number of variables affected (var/node * size) */
1082:   int           *offsets;    /* The canonical variable number for the first variable on each node */
1083:   int           *rows;       /* Rows corresponding to boundary nodes */
1084:   PetscScalar    elem = diag;
1085:   IS             is;
1086:   int            rank;
1087:   int            b, i, j, count;
1088: #ifdef PETSC_USE_BOPT_g
1089:   PetscTruth     opt;
1090: #endif
1091:   int            ierr;

1096:   GMatGetGrid(A, &grid);
1097:   offsets = order->offsets;
1098:   MPI_Comm_rank(grid->comm, &rank);

1100:   /* Allocate memory */
1101:   for(b = 0, totSize = 0, vars = 0; b < num; b++) {
1102:     GridValidField(grid, field[b]);
1103:     GridGetBoundarySize(grid, bd[b], field[b], &size);
1104:     totSize += size;
1105:     vars    += size*grid->fields[field[b]].disc->comp;
1106:   }
1107:   if (totSize == 0) {
1108: #ifdef PETSC_USE_BOPT_g
1109:     PetscSynchronizedFlush(grid->comm);
1110: #endif
1111:     ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is);
1112:     MatZeroRows(A, is, &elem);
1113:     ISDestroy(is);
1114:     return(0);
1115:   }
1116:   PetscMalloc(vars * sizeof(int), &rows);

1118:   /* Loop over boundaries */
1119:   for(b = 0, count = 0; b < num; b++) {
1120:     comp       = grid->fields[field[b]].disc->comp;
1121:     localStart = order->localStart[field[b]];
1122:     /* Loop over boundary nodes */
1123:     GridGetBoundaryStart(grid, bd[b], field[b], PETSC_FALSE, &node, &nclass);
1124:     for(i = 0; node >= 0; i++) {
1125:       for(j = 0; j < comp; j++, count++) {
1126:         rows[count] = offsets[node] + j + localStart[nclass];
1127: #ifdef PETSC_USE_BOPT_g
1128:         PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1129:         if (opt == PETSC_TRUE) {
1130:           PetscSynchronizedPrintf(grid->comm, "[%d]bd %d field: %d node: %d row: %d class: %dn",
1131:                                   rank, bd[b], field[b], node, rows[count], nclass);
1132:         }
1133: #endif
1134:       }
1135:       GridGetBoundaryNext(grid, bd[b], field[b], PETSC_FALSE, &node, &nclass);
1136:     }
1137:   }
1138: #ifdef PETSC_USE_BOPT_g
1139:   PetscSynchronizedFlush(grid->comm);
1140:   if (count != vars) SETERRQ2(PETSC_ERR_PLIB, "Boundary size %d should be %d", count, vars);
1141:   PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
1142: #endif
1143:   /* Set rows of A to the identity */
1144:   ISCreateGeneral(PETSC_COMM_SELF, vars, rows, &is);
1145:   MatZeroRows(A, is, &elem);
1146:   ISDestroy(is);

1148:   GridResetConstrainedMultiply_Private(grid, A);

1150:   PetscFree(rows);
1151:   return(0);
1152: }

1154: #undef  __FUNCT__
1156: /*@C GridSetMatPointBoundary
1157:   This function sets Dirchlet boundary conditions on the linear system matrix arising
1158:   from the underlying grid.

1160:   Collective on GMat

1162:   Input Parameters:
1163: + node  - The constrained node
1164: . field - The field to which the conditions apply
1165: . diag  - The scalar to be placed on the diagonal
1166: - ctx   - The user-supplied context

1168:   Output Parameter:
1169: . A     - The system matrix

1171:   Level: advanced

1173: .keywords boundary conditions, finite element
1174: .seealso MeshGetBoundaryStart
1175: @*/
1176: int GridSetMatPointBoundary(int node, int field, PetscScalar diag, GMat A, void *ctx)
1177: {
1178:   Grid grid;
1179:   int  ierr;

1183:   GMatGetGrid(A, &grid);
1184:   GridSetMatPointBoundaryRectangular(node, field, diag, grid->order, A, ctx);
1185:   return(0);
1186: }

1188: #undef  __FUNCT__
1190: /*@C GridSetMatPointBoundaryRectangular
1191:   This function sets Dirchlet boundary conditions on the linear system matrix arising
1192:   from the underlying grid, and the default variable ordering can be overridden.

1194:   Collective on GMat

1196:   Input Parameters:
1197: + node     - The constrained node
1198: . field    - The field to which the conditions apply
1199: . diag     - The scalar to be placed on the diagonal
1200: . order    - The test variable ordering
1201: - ctx      - The user-supplied context

1203:   Output Parameter:
1204: . A        - The system matrix

1206:   Level: advanced

1208: .keywords boundary conditions, finite element
1209: .seealso MeshGetBoundaryStart
1210: @*/
1211: int GridSetMatPointBoundaryRectangular(int node, int field, PetscScalar diag, VarOrdering order, GMat A, void *ctx)
1212: {
1213:   Grid           grid;
1214:   int            comp;       /* The number of field components */
1215:   int           *localStart; /* The offset of this field on a node of a given class */
1216:   int            nclass;     /* The class of the current boundary node */
1217:   int           *offsets;    /* The canonical variable number for the first variable on each node */
1218:   int           *rows;       /* Rows corresponding to boundary nodes */
1219:   PetscScalar    elem = diag;
1220:   IS             is;
1221:   int            rank;
1222:   int            j;
1223: #ifdef PETSC_USE_BOPT_g
1224:   PetscTruth     opt;
1225: #endif
1226:   int            ierr;

1229:   if (node < 0) {
1230:     ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is);
1231:     MatZeroRows(A, is, &elem);
1232:     ISDestroy(is);
1233:     return(0);
1234:   }
1237:   GMatGetGrid(A, &grid);
1238:   GridValidField(grid, field);
1239:   MPI_Comm_rank(grid->comm, &rank);
1240:   comp       = grid->fields[field].disc->comp;
1241:   offsets    = order->offsets;
1242:   localStart = order->localStart[field];

1244:   /* Allocate memory */
1245:   PetscMalloc(comp * sizeof(int), &rows);

1247:   GridGetNodeClass(grid, node, &nclass);
1248:   for(j = 0; j < comp; j++) {
1249:     rows[j] = offsets[node] + j + localStart[nclass];
1250: #ifdef PETSC_USE_BOPT_g
1251:     PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1252:     if (opt == PETSC_TRUE) {
1253:       PetscPrintf(PETSC_COMM_SELF, "[%d]field: %d node: %d row: %d class: %dn", rank, field, node, rows[j], nclass);
1254:     }
1255: #endif
1256:   }
1257: #ifdef PETSC_USE_BOPT_g
1258:   PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
1259: #endif
1260:   /* Set rows of A to the identity */
1261:   ISCreateGeneral(PETSC_COMM_SELF, comp, rows, &is);
1262:   MatZeroRows(A, is, &elem);
1263:   ISDestroy(is);

1265:   GridResetConstrainedMultiply_Private(grid, A);

1267:   PetscFree(rows);
1268:   return(0);
1269: }

1271: /*------------------------------------------------- Vector Functions ------------------------------------------------*/
1272: #undef  __FUNCT__
1274: /*@C GridSetVecBoundary
1275:   This function sets Dirchlet boundary conditions on the linear Rhs arising
1276:   from the underlying grid.

1278:   Collective on GVec

1280:   Input Parameters:
1281: + bd       - The marker for the boundary to apply conditions along
1282: . field    - The field to which the conditions apply
1283: . f        - The function which defines the boundary condition
1284: - ctx      - The user-supplied context

1286:   Output Parameter:
1287: . b        - The Rhs vector

1289:   Level: advanced

1291: .keywords boundary conditions, finite element
1292: .seealso MeshGetBoundaryStart
1293: @*/
1294: int GridSetVecBoundary(int bd, int field, PointFunction f, GVec b, void *ctx)
1295: {
1296:   Grid grid;
1297:   int  ierr;

1301:   GVecGetGrid(b, &grid);
1302:   GridSetVecBoundaryRectangular(1, &bd, &field, &f, grid->order, b, ctx);
1303:   return(0);
1304: }

1306: #undef  __FUNCT__
1308: /*@C GridSetVecBoundaryRectangular
1309:   This function sets Dirchlet boundary conditions on the linear Rhs arising
1310:   from the underlying grid, and the default variable ordering can be overridden.

1312:   Collective on GVec

1314:   Input Parameters:
1315: + num   - The number of boundary conditions
1316: . bd    - The markers for each boundary to apply conditions along
1317: . field - The fields to which the conditions apply
1318: . f     - The functions which define the boundary conditions
1319: . order - The test variable ordering
1320: - ctx   - The user-supplied context

1322:   Output Parameter:
1323: . b     - The Rhs vector

1325:   Level: advanced

1327: .keywords boundary conditions, finite element
1328: .seealso MeshGetBoundaryStart
1329: @*/
1330: int GridSetVecBoundaryRectangular(int num, int *bd, int *field, PointFunction *f, VarOrdering order, GVec b, void *ctx)
1331: {
1332:   Grid           grid;
1333:   Mesh           mesh;
1334:   int            comp;       /* The number of field components */
1335:   int           *sizes;      /* The number of nodes in each boundary */
1336:   int            totSize;    /* The number of nodes in all boundaries */
1337:   int            maxSize;    /* The maximum number of nodes in any boundary */
1338:   int           *localStart; /* The offset of this field on a node of a given class */
1339:   int            node;       /* The canonical node number of the current boundary node */
1340:   int            nclass;     /* The class of the current boundary node */
1341:   double        *x, *y, *z;  /* Coordinates of the boundary nodes */
1342:   int            vars;       /* The number of variables affected (var/node * size) */
1343:   int           *offsets;    /* The canonical variable number for the first variable on each node */
1344:   int           *rows;       /* Rows corresponding to boundary nodes */
1345:   PetscScalar   *values;     /* Boundary values */
1346:   int            size, rank;
1347:   int            c, i, j, count, off;
1348: #ifdef PETSC_USE_BOPT_g
1349:   PetscTruth     opt;
1350: #endif
1351:   int            ierr;

1356:   GVecGetGrid(b, &grid);
1357:   mesh    = grid->mesh;
1358:   offsets = order->offsets;
1359:   MPI_Comm_rank(grid->comm, &rank);

1361:   /* Support for constrained problems */
1362:   VecGetSize(b, &size);
1363:   if (grid->isConstrained) {
1364:     if (size != grid->constraintOrder->numVars) {
1365:       SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->constraintOrder->numVars);
1366:     }
1367:     offsets = grid->constraintOrder->offsets;
1368:   } else {
1369:     if (size != grid->order->numVars) {
1370:       SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->order->numVars);
1371:     }
1372:   }

1374:   /* Allocate memory */
1375:   PetscMalloc(num * sizeof(int), &sizes);
1376:   for(c = 0, totSize = 0, maxSize = 0, vars = 0; c < num; c++) {
1377:     GridValidField(grid, field[c]);
1378:     GridGetBoundarySize(grid, bd[c], field[c], &sizes[c]);
1379:     totSize += sizes[c];
1380:     maxSize  = PetscMax(maxSize, sizes[c]);
1381:     vars    += sizes[c]*grid->fields[field[c]].disc->comp;
1382:   }
1383:   if (totSize == 0) {
1384: #ifdef PETSC_USE_BOPT_g
1385:     PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1386:     if (opt == PETSC_TRUE) {
1387:       PetscSynchronizedFlush(grid->comm);
1388:       PetscSynchronizedFlush(grid->comm);
1389:     }
1390: #endif
1391:     VecAssemblyBegin(b);
1392:     VecAssemblyEnd(b);
1393:     return(0);
1394:   }
1395:   PetscMalloc(maxSize * sizeof(double),      &x);
1396:   PetscMalloc(maxSize * sizeof(double),      &y);
1397:   PetscMalloc(maxSize * sizeof(double),      &z);
1398:   PetscMalloc(vars    * sizeof(PetscScalar), &values);
1399:   PetscMalloc(vars    * sizeof(int),         &rows);

1401:   /* Loop over boundaries */
1402:   for(c = 0, count = 0, off = 0; c < num; c++, off = count) {
1403:     if (sizes[c] == 0) {
1404: #ifdef PETSC_USE_BOPT_g
1405:     PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1406:     if (opt == PETSC_TRUE) {
1407:         PetscSynchronizedFlush(grid->comm);
1408:         PetscSynchronizedFlush(grid->comm);
1409:     }
1410: #endif
1411:       continue;
1412:     }
1413:     comp       = grid->fields[field[c]].disc->comp;
1414:     localStart = order->localStart[field[c]];
1415:     /* Loop over boundary nodes */
1416:     GridGetBoundaryStart(grid, bd[c], field[c], PETSC_FALSE, &node, &nclass);
1417:     for(i = 0; node >= 0; i++) {
1418:       for(j = 0; j < comp; j++, count++) {
1419:         rows[count] = offsets[node] + j + localStart[nclass];
1420: #ifdef PETSC_USE_BOPT_g
1421:         PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1422:         if (opt == PETSC_TRUE) {
1423:           PetscSynchronizedPrintf(grid->comm, "[%d]bd %d field: %d node: %d row: %d class: %dn",
1424:                                   rank, bd[c], field[c], node, rows[count], nclass);
1425:         }
1426: #endif
1427:       }
1428:       MeshGetNodeCoords(mesh, node, &x[i], &y[i], &z[i]);
1429:       GridGetBoundaryNext(grid, bd[c], field[c], PETSC_FALSE, &node, &nclass);
1430:     }
1431: #ifdef PETSC_USE_BOPT_g
1432:     PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1433:     if (opt == PETSC_TRUE) {
1434:       PetscSynchronizedFlush(grid->comm);
1435:     }
1436: #endif
1437:     /* Get boundary values */
1438:     (*(f[c]))(sizes[c], comp, x, y, z, &values[off], ctx);
1439: #ifdef PETSC_USE_BOPT_g
1440:     PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1441:     if (opt == PETSC_TRUE) {
1442:       PetscPrintf(grid->comm, "Setting boundary values in rhs bd %d field %dn", bd[c], field[c]);
1443:       for(i = off; i < count; i++) PetscSynchronizedPrintf(grid->comm, "  row: %d val: %gn", rows[i], PetscRealPart(values[i]));
1444:       PetscSynchronizedFlush(grid->comm);
1445:     }
1446:     PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
1447: #endif
1448:   }
1449:   if (count != vars) SETERRQ2(PETSC_ERR_PLIB, "Boundary size %d should be %d", count, vars);
1450:   /* Put values in Rhs */
1451:   VecSetValues(b, vars, rows, values, INSERT_VALUES);
1452:   VecAssemblyBegin(b);
1453:   VecAssemblyEnd(b);

1455:   PetscFree(sizes);
1456:   PetscFree(x);
1457:   PetscFree(y);
1458:   PetscFree(z);
1459:   PetscFree(values);
1460:   PetscFree(rows);
1461:   return(0);
1462: }

1464: #undef  __FUNCT__
1466: /*@C GridSetVecPointBoundary
1467:   This function sets Dirchlet boundary conditions on the linear Rhs arising
1468:   from the underlying grid.

1470:   Collective on GVec

1472:   Input Parameters:
1473: + node  - The constrained node
1474: . field - The field to which the conditions apply
1475: . f     - The function which defines the boundary condition
1476: - ctx   - The user-supplied context

1478:   Output Parameter:
1479: . b     - The Rhs vector

1481:   Level: advanced

1483: .keywords boundary conditions, finite element
1484: .seealso MeshGetBoundaryStart
1485: @*/
1486: int GridSetVecPointBoundary(int node, int field, PointFunction f, GVec b, void *ctx)
1487: {
1488:   Grid grid;
1489:   int  ierr;

1493:   GVecGetGrid(b, &grid);
1494:   GridSetVecPointBoundaryRectangular(node, field, f, grid->order, b, ctx);
1495:   return(0);
1496: }

1498: #undef  __FUNCT__
1500: /*@C GridSetVecPointBoundaryRectangular
1501:   This function sets Dirchlet boundary conditions on the linear Rhs arising
1502:   from the underlying grid, and the default variable ordering can be overridden.

1504:   Collective on GVec

1506:   Input Parameters:
1507: + node  - The constriained node
1508: . field - The field to which the conditions apply
1509: . f     - The function which defines the boundary condition
1510: . order - The test variable ordering
1511: - ctx   - The user-supplied context

1513:   Output Parameter:
1514: . b     - The Rhs vector

1516:   Level: advanced

1518: .keywords boundary conditions, finite element
1519: .seealso MeshGetBoundaryStart
1520: @*/
1521: int GridSetVecPointBoundaryRectangular(int node, int field, PointFunction f, VarOrdering order, GVec b, void *ctx)
1522: {
1523:   Grid           grid;
1524:   Mesh           mesh;
1525:   int            comp;       /* The number of field components */
1526:   int            size;       /* The number of nodes in the boundary */
1527:   int           *localStart; /* The offset of this field on a node of a given class */
1528:   int            nclass;     /* The class of the current boundary node */
1529:   double         x, y, z;    /* Coordinates of the boundary nodes */
1530:   int           *offsets;    /* The canonical variable number for the first variable on each node */
1531:   int           *rows;       /* Rows corresponding to boundary nodes */
1532:   PetscScalar   *values;     /* Boundary values */
1533:   int            rank;
1534:   int            c;
1535: #ifdef PETSC_USE_BOPT_g
1536:   PetscTruth     opt;
1537: #endif
1538:   int            ierr;

1541:   if (node < 0) {
1542:     VecAssemblyBegin(b);
1543:     VecAssemblyEnd(b);
1544:     return(0);
1545:   }
1548:   GVecGetGrid(b, &grid);
1549:   GridValidField(grid, field);
1550:   MPI_Comm_rank(grid->comm, &rank);
1551:   mesh       = grid->mesh;
1552:   comp       = grid->fields[field].disc->comp;
1553:   offsets    = order->offsets;
1554:   localStart = order->localStart[field];

1557:   /* Support for constrained problems */
1558:   VecGetSize(b, &size);
1559:   if (grid->isConstrained) {
1560:     if (size != grid->constraintOrder->numVars) {
1561:       SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->constraintOrder->numVars);
1562:     }
1563:     offsets = grid->constraintOrder->offsets;
1564:   } else {
1565:     if (size != grid->order->numVars) {
1566:       SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->order->numVars);
1567:     }
1568:   }

1570:   /* Allocate memory */
1571:   size = 1;
1572:   PetscMalloc(comp * sizeof(PetscScalar), &values);
1573:   PetscMalloc(comp * sizeof(int),         &rows);

1575:   MeshGetNodeCoords(mesh, node, &x, &y, &z);
1576:   GridGetNodeClass(grid, node, &nclass);
1577:   for(c = 0; c < comp; c++) {
1578:     rows[c] = offsets[node] + c + localStart[nclass];
1579: #ifdef PETSC_USE_BOPT_g
1580:     PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1581:     if (opt == PETSC_TRUE) {
1582:       PetscPrintf(PETSC_COMM_SELF, "[%d]field: %d node: %d row: %d class: %dn", rank, field, node, rows[c], nclass);
1583:     }
1584: #endif
1585:   }
1586:   /* Get boundary values */
1587:   (*f)(size, comp, &x, &y, &z, values, ctx);
1588:   /* Put values in Rhs */
1589: #ifdef PETSC_USE_BOPT_g
1590:   PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1591:   if (opt == PETSC_TRUE) {
1592:     PetscPrintf(PETSC_COMM_SELF, "Setting boundary values on rhs node %d field %dn", node, field);
1593:     for(c = 0; c < comp; c++) PetscPrintf(PETSC_COMM_SELF, "  row: %d val: %gn", rows[c], PetscRealPart(values[c]));
1594:   }
1595:   PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
1596: #endif
1597:   VecSetValues(b, comp, rows, values, INSERT_VALUES);
1598:   VecAssemblyBegin(b);
1599:   VecAssemblyEnd(b);

1601:   PetscFree(values);
1602:   PetscFree(rows);
1603:   return(0);
1604: }

1606: #undef  __FUNCT__
1608: /*@C GridSetVecBoundaryDifference
1609:   This function sets Dirchlet boundary conditions on the linear Rhs arising
1610:   from the underlying grid, but actually sets it to the difference of the
1611:   function value and the value in the given vector. This is commonly used in
1612:   a time dependent, nonlinear problem for which we would like the rhs boundary
1613:   values to be:

1615:       U^{n+1}_k - U^{n+1}_{k+1}

1617:   where n is the time iteration index, and k is the Newton iteration index. This
1618:   means that the solution will be updated to U^{n+1}_{k+1} if the Jacobian is the
1619:   identity for that row. This is very useful for time dependent boundary conditions
1620:   for which the traditional method of letting the rhs value be zero does not work.

1622:   Collective on GVec

1624:   Input Parameters:
1625: + bd       - The marker for the boundary to apply conditions along
1626: . u        - A grid vector, usually the previous solution
1627: . field    - The field to which the conditions apply
1628: . f        - The function which defines the boundary condition
1629: - ctx      - The user-supplied context

1631:   Output Parameter:
1632: . b        - The Rhs vector

1634:   Level: advanced

1636: .keywords boundary conditions, finite element
1637: .seealso MeshGetBoundaryStart
1638: @*/
1639: int GridSetVecBoundaryDifference(int bd, int field, GVec u, PointFunction f, GVec b, void *ctx)
1640: {
1641:   Grid grid;
1642:   int  ierr;

1646:   GVecGetGrid(b, &grid);
1647:   GridSetVecBoundaryDifferenceRectangular(bd, field, u, f, grid->order, b, ctx);
1648:   return(0);
1649: }

1651: #undef  __FUNCT__
1653: /*@C GridSetVecBoundaryDifferenceRectangular
1654:   This function sets Dirchlet boundary conditions on the linear Rhs arising
1655:   from the underlying grid, but actually sets it to the difference of the
1656:   function value and the value in the given vector. This is commonly used in
1657:   a time dependent, nonlinear problem for which we would like the rhs boundary
1658:   values to be:

1660:       U^{n+1}_k - U^{n+1}_{k+1}

1662:   where n is the time iteration index, and k is the Newton iteration index. This
1663:   means that the solution will be updated to U^{n+1}_{k+1} if the Jacobian is the
1664:   identity for that row. This is very useful for time dependent boundary conditions
1665:   for which the traditional method of letting the rhs value be zero does not work.

1667:   Collective on GVec

1669:   Input Parameters:
1670: + bd    - The marker for the boundary to apply conditions along
1671: . u     - A grid vector, usually the previous solution
1672: . field - The field to which the conditions apply
1673: . f     - The function which defines the boundary condition
1674: . order - The test variable ordering
1675: - ctx   - The user-supplied context

1677:   Output Parameter:
1678: . b     - The Rhs vector

1680:   Level: advanced

1682: .keywords boundary conditions, finite element
1683: .seealso MeshGetBoundaryStart
1684: @*/
1685: int GridSetVecBoundaryDifferenceRectangular(int bd, int field, GVec u, PointFunction f, VarOrdering order, GVec b, void *ctx)
1686: {
1687:   Grid           grid;
1688:   Mesh           mesh;
1689:   int            comp;       /* The number of field components */
1690:   int            size;       /* The number of nodes in the boundary */
1691:   int           *localStart; /* The offset of this field on a node of a given class */
1692:   int            node;       /* The canonical node number of the current boundary node */
1693:   int            nclass;     /* The class of the current boundary node */
1694:   double        *x, *y, *z;  /* Coordinates of the boundary nodes */
1695:   int            vars;       /* The number of variables affected (var/node * size) */
1696:   int           *offsets;    /* The canonical variable number for the first variable on each node */
1697:   int           *rows;       /* Rows corresponding to boundary nodes */
1698:   PetscScalar   *values;     /* Boundary values */
1699:   PetscScalar   *uArray;     /* The values in the vector u */
1700:   int            firstVar;   /* The canonical number of the first variable in this domain */
1701:   int            rank;
1702:   int            i, j, count;
1703: #ifdef PETSC_USE_BOPT_g
1704:   PetscTruth     opt;
1705: #endif
1706:   int            ierr;

1711:   GVecGetGrid(b, &grid);
1712:   GridValidField(grid, field);
1713:   MPI_Comm_rank(grid->comm, &rank);
1714:   mesh       = grid->mesh;
1715:   comp       = grid->fields[field].disc->comp;
1716:   firstVar   = order->firstVar[rank];
1717:   offsets    = order->offsets;
1718:   localStart = order->localStart[field];

1720:   /* Support for constrained problems */
1721:   VecGetSize(b, &size);
1722:   if (grid->isConstrained) {
1723:     if (size != grid->constraintOrder->numVars) {
1724:       SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->constraintOrder->numVars);
1725:     }
1726:     offsets = grid->constraintOrder->offsets;
1727:   } else {
1728:     if (size != grid->order->numVars) {
1729:       SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->order->numVars);
1730:     }
1731:   }

1733:   /* Allocate memory */
1734:   GridGetBoundarySize(grid, bd, field, &size);
1735:   if (size == 0) {
1736: #ifdef PETSC_USE_BOPT_g
1737:     PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1738:     if (opt == PETSC_TRUE) {
1739:       PetscSynchronizedFlush(grid->comm);
1740:       PetscSynchronizedFlush(grid->comm);
1741:     }
1742: #endif
1743:     VecAssemblyBegin(b);
1744:     VecAssemblyEnd(b);
1745:     return(0);
1746:   }
1747:   vars = size*comp;
1748:   PetscMalloc(size * sizeof(double),      &x);
1749:   PetscMalloc(size * sizeof(double),      &y);
1750:   PetscMalloc(size * sizeof(double),      &z);
1751:   PetscMalloc(vars * sizeof(PetscScalar), &values);
1752:   PetscMalloc(vars * sizeof(int),         &rows);

1754:   /* Loop over boundary nodes */
1755:   GridGetBoundaryStart(grid, bd, field, PETSC_FALSE, &node, &nclass);
1756:   for(i = 0, count = 0; node >= 0; i++) {
1757:     for(j = 0; j < comp; j++, count++) {
1758:       rows[count] = offsets[node] + j + localStart[nclass];
1759: #ifdef PETSC_USE_BOPT_g
1760:       PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1761:       if (opt == PETSC_TRUE) {
1762:         PetscSynchronizedPrintf(grid->comm, "[%d]bd %d field: %d node: %d row: %d class: %dn",
1763:                                 rank, bd, field, node, rows[count], nclass);
1764:       }
1765: #endif
1766:     }
1767:     MeshGetNodeCoords(mesh, node, &x[i], &y[i], &z[i]);
1768:     GridGetBoundaryNext(grid, bd, field, PETSC_FALSE, &node, &nclass);
1769:   }
1770: #ifdef PETSC_USE_BOPT_g
1771:   PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1772:   if (opt == PETSC_TRUE) {
1773:     PetscSynchronizedFlush(grid->comm);
1774:   }
1775: #endif
1776:   /* Get boundary values */
1777:   (*f)(size, comp, x, y, z, values, ctx);
1778:   /* Taking the difference (we know that no values are off-processor) */
1779:   VecGetArray(u, &uArray);
1780:   for(i = 0; i < vars; i++)
1781:     values[i] =  uArray[rows[i]-firstVar] - values[i];
1782:   VecRestoreArray(u, &uArray);
1783:   /* Put values in Rhs */
1784: #ifdef PETSC_USE_BOPT_g
1785:   PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1786:   if (opt == PETSC_TRUE) {
1787:     PetscPrintf(grid->comm, "Setting boundary values in rhs bd %d field %dn", bd, field);
1788:     for(i = 0; i < vars; i++) PetscSynchronizedPrintf(grid->comm, "  row: %d val: %gn", rows[i], PetscRealPart(values[i]));
1789:     PetscSynchronizedFlush(grid->comm);
1790:   }
1791:   PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
1792: #endif
1793:   VecSetValues(b, vars, rows, values, INSERT_VALUES);
1794:   VecAssemblyBegin(b);
1795:   VecAssemblyEnd(b);

1797:   PetscFree(x);
1798:   PetscFree(y);
1799:   PetscFree(z);
1800:   PetscFree(values);
1801:   PetscFree(rows);
1802:   return(0);
1803: }

1805: #undef  __FUNCT__
1807: /*@C GridSetVecPointBoundaryDifference
1808:   This function sets Dirchlet boundary conditions on the linear Rhs arising
1809:   from the underlying grid, but actually sets it to the difference of the
1810:   function value and the value in the given vector. This is commonly used in
1811:   a time dependent, nonlinear problem for which we would like the rhs boundary
1812:   values to be:

1814:       U^{n+1}_k - U^{n+1}_{k+1}

1816:   where n is the time iteration index, and k is the Newton iteration index. This
1817:   means that the solution will be updated to U^{n+1}_{k+1} if the Jacobian is the
1818:   identity for that row. This is very useful for time dependent boundary conditions
1819:   for which the traditional method of letting the rhs value be zero does not work.

1821:   Collective on GVec

1823:   Input Parameters:
1824: + node     - The constrained node
1825: . u        - A grid vector, usually the previous solution
1826: . field    - The field to which the conditions apply
1827: . f        - The function which defines the boundary condition
1828: - ctx      - The user-supplied context

1830:   Output Parameter:
1831: . b        - The Rhs vector

1833:   Level: advanced

1835: .keywords boundary conditions, finite element
1836: .seealso MeshGetBoundaryStart
1837: @*/
1838: int GridSetVecPointBoundaryDifference(int node, int field, GVec u, PointFunction f, GVec b, void *ctx)
1839: {
1840:   Grid grid;
1841:   int  ierr;

1845:   GVecGetGrid(b, &grid);
1846:   GridSetVecPointBoundaryDifferenceRectangular(node, field, u, f, grid->order, b, ctx);
1847:   return(0);
1848: }

1850: #undef  __FUNCT__
1852: /*@C GridSetVecBoundaryDifferenceRectangular
1853:   This function sets Dirchlet boundary conditions on the linear Rhs arising
1854:   from the underlying grid, but actually sets it to the difference of the
1855:   function value and the value in the given vector. This is commonly used in
1856:   a time dependent, nonlinear problem for which we would like the rhs boundary
1857:   values to be:

1859:       U^{n+1}_k - U^{n+1}_{k+1}

1861:   where n is the time iteration index, and k is the Newton iteration index. This
1862:   means that the solution will be updated to U^{n+1}_{k+1} if the Jacobian is the
1863:   identity for that row. This is very useful for time dependent boundary conditions
1864:   for which the traditional method of letting the rhs value be zero does not work.

1866:   Collective on GVec

1868:   Input Parameters:
1869: + node  - The constrained node
1870: . u     - A grid vector, usually the previous solution
1871: . field - The field to which the conditions apply
1872: . f     - The function which defines the boundary condition
1873: . order - The test variable ordering
1874: - ctx   - The user-supplied context

1876:   Output Parameter:
1877: . b     - The Rhs vector

1879:   Level: advanced

1881: .keywords boundary conditions, finite element
1882: .seealso MeshGetBoundaryStart
1883: @*/
1884: int GridSetVecPointBoundaryDifferenceRectangular(int node, int field, GVec u, PointFunction f, VarOrdering order, GVec b, void *ctx)
1885: {
1886:   Grid           grid;
1887:   Mesh           mesh;
1888:   int            comp;       /* The number of field components */
1889:   int            size;       /* The number of nodes in the boundary */
1890:   int           *localStart; /* The offset of this field on a node of a given class */
1891:   int            nclass;     /* The class of the current boundary node */
1892:   double         x, y, z;    /* Coordinates of the boundary nodes */
1893:   int           *offsets;    /* The canonical variable number for the first variable on each node */
1894:   int           *rows;       /* Rows corresponding to boundary nodes */
1895:   PetscScalar   *values;     /* Boundary values */
1896:   PetscScalar   *uArray;     /* The values in the vector u */
1897:   int            firstVar;   /* The canonical number of the first variable in this domain */
1898:   int            rank;
1899:   int            i, j;
1900: #ifdef PETSC_USE_BOPT_g
1901:   PetscTruth     opt;
1902: #endif
1903:   int            ierr;

1906:   if (node < 0) {
1907:     VecAssemblyBegin(b);
1908:     VecAssemblyEnd(b);
1909:     return(0);
1910:   }
1912:   GVecGetGrid(b, &grid);
1913:   GridGetMesh(grid, &mesh);
1914:   GridValidField(grid, field);
1915:   MPI_Comm_rank(grid->comm, &rank);
1916:   comp       = grid->fields[field].disc->comp;
1917:   firstVar   = order->firstVar[rank];
1918:   offsets    = order->offsets;
1919:   localStart = order->localStart[field];

1921:   /* Support for constrained problems */
1922:   VecGetSize(b, &size);
1923:   if (grid->isConstrained) {
1924:     if (size != grid->constraintOrder->numVars) {
1925:       SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->constraintOrder->numVars);
1926:     }
1927:     offsets = grid->constraintOrder->offsets;
1928:   } else {
1929:     if (size != grid->order->numVars) {
1930:       SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid vector size %d should be %d", size, grid->order->numVars);
1931:     }
1932:   }

1934:   /* Allocate memory */
1935:   size = 1;
1936:   PetscMalloc(comp * sizeof(PetscScalar), &values);
1937:   PetscMalloc(comp * sizeof(int),         &rows);

1939:   MeshGetNodeCoords(mesh, node, &x, &y, &z);
1940:   GridGetNodeClass(grid, node, &nclass);
1941:   for(j = 0; j < comp; j++) {
1942:     rows[j] = offsets[node] + j + localStart[nclass];
1943: #ifdef PETSC_USE_BOPT_g
1944:     PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1945:     if (opt == PETSC_TRUE) {
1946:       PetscPrintf(PETSC_COMM_SELF, "[%d]field: %d node: %d row: %d class: %dn", rank, field, node, rows[j], nclass);
1947:     }
1948: #endif
1949:   }
1950:   /* Get boundary values */
1951:   (*f)(size, comp, &x, &y, &z, values, ctx);
1952:   /* Taking the difference (we know that no values are off-processor) */
1953:   VecGetArray(u, &uArray);
1954:   for(i = 0; i < comp; i++) values[i] =  uArray[rows[i]-firstVar] - values[i];
1955:   VecRestoreArray(u, &uArray);
1956:   /* Put values in Rhs */
1957: #ifdef PETSC_USE_BOPT_g
1958:   PetscOptionsHasName(PETSC_NULL, "-trace_bc", &opt);
1959:   if (opt == PETSC_TRUE) {
1960:     PetscPrintf(grid->comm, "Setting boundary values on rhs node %d field %dn", node, field);
1961:     for(i = 0; i < comp; i++) PetscPrintf(PETSC_COMM_SELF, "  row: %d val: %gn", rows[i], PetscRealPart(values[i]));
1962:   }
1963:   PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
1964: #endif
1965:   VecSetValues(b, comp, rows, values, INSERT_VALUES);
1966:   VecAssemblyBegin(b);
1967:   VecAssemblyEnd(b);

1969:   PetscFree(values);
1970:   PetscFree(rows);
1971:   return(0);
1972: }