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: }