Actual source code: gridDB.c
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: gridDB.c,v 1.4 2000/07/16 05:40:27 knepley Exp $";
3: #endif
5: #include src/grid/gridimpl.h
7: /*-------------------------------------------------- Field Functions ------------------------------------------------*/
8: #undef __FUNCT__
10: /*@C
11: GridAddField - This function defines a field on the grid.
13: Not collective
15: Input Parameters:
16: + grid - The grid
17: . name - The field name
18: . discType - The DiscretizationType for the field
19: - numComp - The nubmer of components in the field
21: Output Parameter:
22: . field - The canonical field number
24: Level: intermediate
26: .keywords: grid, field,
27: .seealso: GridSetFieldName()
28: @*/
29: int GridAddField(Grid grid, const char name[], DiscretizationType discType, int numComp, int *field)
30: {
31: Field *tempFields;
32: int ierr;
36: if (numComp <= 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of field components %d", numComp);
37: while (grid->numFields >= grid->maxFields) {
38: PetscMalloc(grid->maxFields*2 * sizeof(Field), &tempFields);
39: PetscLogObjectMemory(grid, grid->maxFields * sizeof(Field));
40: PetscMemcpy(tempFields, grid->fields, grid->maxFields * sizeof(Field));
41: PetscFree(grid->fields);
42: grid->fields = tempFields;
43: grid->maxFields *= 2;
44: }
45: if (field != PETSC_NULL) {
47: *field = grid->numFields;
48: }
50: PetscStrallocpy(name, &grid->fields[grid->numFields].name);
51: PetscStrallocpy(discType, &grid->fields[grid->numFields].discType);
52: grid->fields[grid->numFields].numComp = numComp;
53: grid->fields[grid->numFields].isConstrained = PETSC_FALSE;
54: grid->fields[grid->numFields].constraintCompDiff = 0;
55: grid->fields[grid->numFields].isActive = PETSC_FALSE;
56: DiscretizationCreate(grid->comm, &grid->fields[grid->numFields].disc);
57: DiscretizationSetNumComponents(grid->fields[grid->numFields].disc, numComp);
58: DiscretizationSetField(grid->fields[grid->numFields].disc, grid->numFields);
59: DiscretizationSetType(grid->fields[grid->numFields].disc, discType);
60: DiscretizationSetFromOptions(grid->fields[grid->numFields].disc);
61: grid->numFields++;
62: return(0);
63: }
65: #undef __FUNCT__
67: /*@C
68: GridGetNumFields - This function gets the number of fields in the grid.
70: Not collective
72: Input Parameter:
73: . grid - The grid
75: Output Parameter:
76: . num - The number of fields
78: Level: intermediate
80: .keywords: grid, field,
81: .seealso: GridSetFieldName()
82: @*/
83: int GridGetNumFields(Grid grid, int *num)
84: {
88: *num = grid->numFields;
89: return(0);
90: }
92: #undef __FUNCT__
94: /*@C
95: GridGetFieldName - This function gets the name of a particular field.
97: Not collective
99: Input Parameters:
100: + grid - The grid
101: - field - The field to name
103: Output Parameter:
104: . name - The name string
106: Level: intermediate
108: .keywords: grid, field, name
109: .seealso: GridSetFieldName()
110: @*/
111: int GridGetFieldName(Grid grid, int field, char **name)
112: {
118: GridValidField(grid, field);
119: if (grid->fields[field].name == PETSC_NULL) {
120: *name = PETSC_NULL;
121: } else {
122: PetscStrallocpy(grid->fields[field].name, name);
123: }
124: return(0);
125: }
127: #undef __FUNCT__
129: /*@C
130: GridSetFieldName - This function sets the name of a particular field.
132: Collective on Grid
134: Input Parameters:
135: + grid - The grid
136: . field - The field to name
137: - name - The name string
139: Level: intermediate
141: .keywords: grid, field, name
142: .seealso: GridGetFieldName()
143: @*/
144: int GridSetFieldName(Grid grid, int field, char *name)
145: {
150: GridValidField(grid, field);
152: PetscStrfree(grid->fields[field].name);
153: PetscStrallocpy(name, &grid->fields[field].name);
154: return(0);
155: }
157: #undef __FUNCT__
159: /*@C
160: GridGetFieldComponents - This function gets the number of components in a field.
162: Not collective
164: Input Parameters
165: + grid - The grid
166: - field - The field
168: Output Parameter:
169: . comp - The number of components in the field
171: Level: intermediate
173: .keywords: grid, field, component
174: .seealso: GridGetFieldDiscType()
175: @*/
176: int GridGetFieldComponents(Grid grid, int field, int *comp)
177: {
180: GridValidField(grid, field);
182: *comp = grid->fields[field].numComp;
183: return(0);
184: }
186: #undef __FUNCT__
188: /*@C
189: GridGetFieldDisc - This function gets the discretization for a field.
191: Not collective
193: Input Parameters:
194: + grid - The grid
195: - field - The field
197: Output Parameter:
198: . disc - The discretization
200: Level: intermediate
202: .keywords: grid, field, discretization
203: .seealso: GridGetFieldComponents()
204: @*/
205: int GridGetFieldDisc(Grid grid, int field, Discretization *disc)
206: {
210: GridValidField(grid, field);
211: *disc = grid->fields[field].disc;
212: return(0);
213: }
215: #undef __FUNCT__
217: /*@C
218: GridGetNumActiveFields - This function gets the number of fields participating in the calculation.
220: Not collective
222: Input Parameter:
223: . grid - The grid
225: Output Parameter:
226: . num - The number of active fields
228: Level: intermediate
230: .keywords: grid, field, active
231: .seealso: GridSetFieldName()
232: @*/
233: int GridGetNumActiveFields(Grid grid, int *num)
234: {
238: *num = grid->numActiveFields;
239: return(0);
240: }
242: #undef __FUNCT__
244: /*@C
245: GridGetActiveField - This function gets the canonical field number for a certain active field.
247: Not collective
249: Input Parameters:
250: + grid - The grid
251: - n - The nth active field
253: Output Parameter:
254: . field - The field
256: Level: intermediate
258: .keywords: grid, field, active
259: .seealso: GridSetFieldName()
260: @*/
261: int GridGetActiveField(Grid grid, int n, int *field)
262: {
266: if ((n < 0) || (n >= grid->numActiveFields)) {
267: SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid active field number %d should be in [0,%d)", n, grid->numActiveFields);
268: }
269: *field = grid->defaultFields[n];
270: return(0);
271: }
273: #undef __FUNCT__
275: /*@C
276: GridSetActiveField - This function makes one field active for a particular grid.
278: Collective on Grid
280: Input Parameters:
281: + grid - The grid
282: - field - The field to make active
284: Level: intermediate
286: .keywords: grid, active, field
287: .seealso: GridAddActiveField()
288: @*/
289: int GridSetActiveField(Grid grid, int field)
290: {
291: int fieldIndex;
296: GridValidField(grid, field);
297: grid->numActiveFields = 0;
298: for(fieldIndex = 0; fieldIndex < grid->numFields; fieldIndex++) grid->fields[fieldIndex].isActive = PETSC_FALSE;
299: GridAddActiveField(grid, field);
300: return(0);
301: }
303: #undef __FUNCT__
305: /*@C
306: GridAddActiveField - This function makes another field active for a particular grid.
308: Collective on Grid
310: Input Parameters:
311: + grid - The grid
312: - field - The field to make active
314: Level: intermediate
316: .keywords: grid, active, field
317: .seealso: GridSetActiveField()
318: @*/
319: int GridAddActiveField(Grid grid, int field)
320: {
321: int *tempFields;
322: int ierr;
326: GridValidField(grid, field);
327: /* Setup assembly variables */
328: while (grid->numActiveFields >= grid->maxActiveFields) {
329: PetscMalloc(grid->maxActiveFields*2 * sizeof(int), &tempFields);
330: PetscLogObjectMemory(grid, grid->maxActiveFields * sizeof(int));
331: PetscMemcpy(tempFields, grid->defaultFields, grid->maxActiveFields * sizeof(int));
332: PetscFree(grid->defaultFields);
333: grid->defaultFields = tempFields;
334: grid->maxActiveFields *= 2;
335: }
336: if (grid->fields[field].isActive == PETSC_FALSE) grid->defaultFields[grid->numActiveFields++] = field;
337: grid->fields[field].isActive = PETSC_TRUE;
338: grid->setupcalled = PETSC_FALSE;
339: grid->bdSetupCalled = PETSC_FALSE;
340: return(0);
341: }
343: /*------------------------------------------------ Variable Functions -----------------------------------------------*/
344: #undef __FUNCT__
346: /*@C GridGetNumVars
347: This function gets the number of variables in the grid.
349: Not collective
351: Input Parameter:
352: . grid - The grid
354: Output Parameters:
355: + locNum - The number of local variables
356: . num - The number of variables
357: . locConstNum - The number of local constrained variables
358: - constNum - The number of constrained variables
360: Level: intermediate
362: .keywords grid, field,
363: .seealso GridSetFieldName
364: @*/
365: int GridGetNumVars(Grid grid, int *locNum, int *num, int *locConstNum, int *constNum)
366: {
371: if (grid->order == PETSC_NULL) {
372: GridSetUp(grid);
373: }
374: if (locNum != PETSC_NULL) {
376: *locNum = grid->order->numLocVars;
377: }
378: if (num != PETSC_NULL) {
380: *num = grid->order->numVars;
381: }
382: if (grid->isConstrained == PETSC_TRUE) {
383: if ((grid->constraintOrder == PETSC_NULL) && ((locConstNum != PETSC_NULL) || (constNum != PETSC_NULL))) {
384: SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Constraints are not yet setup");
385: }
386: if (locConstNum != PETSC_NULL) {
388: *locConstNum = grid->constraintOrder->numLocVars;
389: }
390: if (constNum != PETSC_NULL) {
392: *constNum = grid->constraintOrder->numVars;
393: }
394: }
395: return(0);
396: }
398: /*-------------------------------------------------- Rhs Functions --------------------------------------------------*/
399: #undef __FUNCT__
401: /*@
402: GridGetNumNonlinearOperators - Returns the number of nonlinear operators for the problem defined on the grid.
404: Not collective
406: Input Parameter:
407: . grid - The grid
409: Output Parameter:
410: . numOps - The number of operators
412: Level: intermediate
414: .keywords: grid, operator, nonlinear
415: .seealso: GridGetNonlinearOperatorStart()
416: @*/
417: int GridGetNumNonlinearOperators(Grid grid, int *numOps)
418: {
419: int op;
424: *numOps = 0;
425: for(op = 0; op < grid->numRhsOps; op++) {
426: if (grid->rhsOps[op].nonlinearOp != PETSC_NULL) *numOps += 1;
427: }
428: return(0);
429: }
431: #undef __FUNCT__
433: /*@
434: GridGetNonlinearOperatorStart - Retrieves information about the first nonlinear
435: operator for the problem defined on the grid.
437: Not collective
439: Input Parameter:
440: . grid - The grid
442: Output Parameters:
443: + op - The operator
444: . field - The shape function field
445: . alpha - The multiplier of this operator
446: - isALE - The flag for ALE operators
448: Level: intermediate
450: .keywords: grid, operator, nonlinear
451: .seealso: GridGetNonlinearOperatorNext(), GridGetNumNonlinearOperators()
452: @*/
453: int GridGetNonlinearOperatorStart(Grid grid, NonlinearOperator *op, int *field, PetscScalar *alpha, PetscTruth *isALE)
454: {
459: grid->activeNonlinearOp = -1;
460: GridGetNonlinearOperatorNext(grid, op, field, alpha, isALE);
461: return(0);
462: }
464: #undef __FUNCT__
466: /*@
467: GridGetNonlinearOperatorNext - Retrieves information about the next nonlinear
468: operator for the problem defined on the grid.
470: Not collective
472: Input Parameter:
473: . grid - The grid
475: Output Parameters:
476: + op - The operator
477: . field - The shape function field
478: . alpha - The multiplier of this operator
479: - isALE - The flag for ALE operators
481: Level: intermediate
483: .keywords: grid, operator, nonlinear
484: .seealso: GridNonlinearGetOperatorStart(), GridGetNumNonlinearOperators()
485: @*/
486: int GridGetNonlinearOperatorNext(Grid grid, NonlinearOperator *op, int *field, PetscScalar *alpha, PetscTruth *isALE)
487: {
490: grid->activeNonlinearOp++;
491: while((grid->rhsOps[grid->activeNonlinearOp].nonlinearOp == PETSC_NULL) && (grid->activeNonlinearOp < grid->numRhsOps)) {
492: grid->activeNonlinearOp++;
493: }
494: if (grid->activeNonlinearOp >= grid->numRhsOps) {
495: if (op != PETSC_NULL) *op = PETSC_NULL;
496: if (field != PETSC_NULL) *field = -1;
497: if (alpha != PETSC_NULL) *alpha = 0.0;
498: if (isALE != PETSC_NULL) *isALE = PETSC_FALSE;
499: grid->activeNonlinearOp = -1;
500: return(0);
501: }
502: if (op != PETSC_NULL) {
504: *op = grid->rhsOps[grid->activeNonlinearOp].nonlinearOp;
505: }
506: if (field != PETSC_NULL) {
508: *field = grid->rhsOps[grid->activeNonlinearOp].field;
509: }
510: if (alpha != PETSC_NULL) {
512: *alpha = grid->rhsOps[grid->activeNonlinearOp].alpha;
513: }
514: if (isALE != PETSC_NULL) {
516: *isALE = grid->rhsOps[grid->activeNonlinearOp].isALE;
517: }
518: return(0);
519: }
521: #undef __FUNCT__
523: /*@C GridAddRhsFunction
524: This function adds the point function to use on the rhs for the given field.
526: Collective on Grid
528: Input Parameters:
529: + grid - The grid
530: . field - The field
531: . func - The point function
532: - alpha - A scalar multiplier
534: Level: beginner
536: .keywords active field
537: .seealso GridSetMatOperator, GridSetRhsOperator
538: @*/
539: int GridAddRhsFunction(Grid grid, int field, PointFunction f, PetscScalar alpha)
540: {
541: GridFunc *tempFuncs;
542: int ierr;
546: GridValidField(grid, field);
547: while (grid->numRhsFuncs >= grid->maxRhsFuncs) {
548: PetscMalloc(grid->maxRhsFuncs*2 * sizeof(GridOp), &tempFuncs);
549: PetscLogObjectMemory(grid, grid->maxRhsFuncs * sizeof(GridOp));
550: PetscMemcpy(tempFuncs, grid->rhsFuncs, grid->maxRhsFuncs * sizeof(GridOp));
551: PetscFree(grid->rhsFuncs);
552: grid->rhsFuncs = tempFuncs;
553: grid->maxRhsFuncs *= 2;
554: }
555: grid->rhsFuncs[grid->numRhsFuncs].func = f;
556: grid->rhsFuncs[grid->numRhsFuncs].field = field;
557: grid->rhsFuncs[grid->numRhsFuncs].alpha = alpha;
558: grid->numRhsFuncs++;
559: return(0);
560: }
562: #undef __FUNCT__
564: /*@C GridAddRhsNonlinearOperator
565: This function add the nonlinear operator to use on the Rhs for the given field.
567: Collective on Grid
569: Input Parameters:
570: + grid - The grid
571: . field - The field
572: . f - The function defining the nonlinear operator
573: . alpha - A scalar multiplier
574: - isALE - A flag for ALE operators
576: Level: beginner
578: .keywords active field
579: .seealso GridSetRhsOperator, GridAddRhsOperator
580: @*/
581: int GridAddRhsNonlinearOperator(Grid grid, int field, NonlinearOperator f, PetscScalar alpha, PetscTruth isALE)
582: {
583: GridOp *tempOps;
584: int ierr;
588: GridValidField(grid, field);
589: while (grid->numRhsOps >= grid->maxRhsOps) {
590: PetscMalloc(grid->maxRhsOps*2 * sizeof(GridOp), &tempOps);
591: PetscLogObjectMemory(grid, grid->maxRhsOps * sizeof(GridOp));
592: PetscMemcpy(tempOps, grid->rhsOps, grid->maxRhsOps * sizeof(GridOp));
593: PetscFree(grid->rhsOps);
594: grid->rhsOps = tempOps;
595: grid->maxRhsOps *= 2;
596: }
597: grid->rhsOps[grid->numRhsOps].nonlinearOp = f;
598: grid->rhsOps[grid->numRhsOps].field = field;
599: grid->rhsOps[grid->numRhsOps].alpha = alpha;
600: grid->rhsOps[grid->numRhsOps].isALE = isALE;
601: grid->rhsOps[grid->numRhsOps].op = -1;
602: grid->numRhsOps++;
603: if (isALE == PETSC_TRUE) grid->ALEActive = PETSC_TRUE;
604: return(0);
605: }
607: #undef __FUNCT__
609: /*@C GridSetRhsOperator
610: This function sets the operator to use on the Rhs for the given field.
612: Collective on Grid
614: Input Parameters:
615: + grid - The grid
616: . sfield - The field for shape functions
617: . tfield - The field for test functions
618: . op - The operator
619: . alpha - A scalar multiple for the operator
620: - isALE - A flag for ALE operators
622: Output Parameter:
623: . index - [Optional] A unique indentifier for the operator
625: Level: beginner
627: .keywords active field
628: .seealso GridAddRhsOperator, GridAddMatOperator, GridSetRhsFunction
629: @*/
630: int GridSetRhsOperator(Grid grid, int sField, int tField, int op, PetscScalar alpha, PetscTruth isALE, int *index)
631: {
636: grid->numRhsOps = 0;
637: GridAddRhsOperator(grid, sField, tField, op, alpha, isALE, index);
638: return(0);
639: }
641: #undef __FUNCT__
643: /*@C GridAddRhsOperator
644: This function adds the operator to the Rhs for the given field.
646: Collective on Grid
648: Input Parameters:
649: + grid - The grid
650: . sfield - The field for shape functions
651: . tfield - The field for test functions
652: . op - The operator
653: . alpha - A scalar multiple for the operator
654: - isALE - A flag for ALE operators
656: Output Parameter:
657: . index - [Optional] A unique indentifier for the operator
659: Level: beginner
661: .keywords active field
662: .seealso GridSetRhsOperator, GridSetMatOperator, GridSetRhsFunction
663: @*/
664: int GridAddRhsOperator(Grid grid, int sField, int tField, int op, PetscScalar alpha, PetscTruth isALE, int *index)
665: {
666: Discretization sDisc = grid->fields[sField].disc;
667: Discretization tDisc = grid->fields[tField].disc;
668: int dim = grid->dim;
669: GridOp *tempOps;
670: int ierr;
674: GridValidField(grid, sField);
675: GridValidField(grid, tField);
676: if ((op < 0) || (op >= sDisc->numOps) || (!sDisc->operators[op])) {
677: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d", op);
678: }
679: if ((sDisc->operators[op]->precompInt == PETSC_NULL) &&
680: (sDisc->operators[op]->opFunc == PETSC_NULL) &&
681: (sDisc->operators[op]->ALEOpFunc == PETSC_NULL)) {
682: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d has no associated function", op);
683: }
684: if (sDisc->numQuadPoints != tDisc->numQuadPoints) {
685: SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible quadrature sets");
686: }
687: while (grid->numRhsOps >= grid->maxRhsOps) {
688: PetscMalloc(grid->maxRhsOps*2 * sizeof(GridOp), &tempOps);
689: PetscLogObjectMemory(grid, grid->maxRhsOps * sizeof(GridOp));
690: PetscMemcpy(tempOps, grid->rhsOps, grid->maxRhsOps * sizeof(GridOp));
691: PetscFree(grid->rhsOps);
692: grid->rhsOps = tempOps;
693: grid->maxRhsOps *= 2;
694: }
695: if (index != PETSC_NULL) {
697: *index = grid->numRhsOps;
698: }
699: grid->rhsOps[grid->numRhsOps].op = op;
700: grid->rhsOps[grid->numRhsOps].field = sField;
701: grid->rhsOps[grid->numRhsOps].alpha = alpha;
702: grid->rhsOps[grid->numRhsOps].isALE = isALE;
703: grid->rhsOps[grid->numRhsOps].nonlinearOp = PETSC_NULL;
704: grid->numRhsOps++;
705: if (isALE == PETSC_TRUE) grid->ALEActive = PETSC_TRUE;
706: /* Assign test function field to operator -- Need more sophisticated checking */
707: switch (op) {
708: case IDENTITY:
709: case LAPLACIAN:
710: if (sDisc->size != tDisc->size) {
711: SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
712: }
713: break;
714: case GRADIENT:
715: if (sDisc->comp*dim != tDisc->comp) {
716: SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
717: }
718: break;
719: case DIVERGENCE:
720: if (tDisc->comp*dim != sDisc->comp) {
721: SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
722: }
723: break;
724: }
725: sDisc->operators[op]->test = tDisc;
726: return(0);
727: }
729: #undef __FUNCT__
731: /*@C GridScaleRhs
732: This function multiplies all the current components of the Rhs by the same scalar.
734: Collective on Grid
736: Input Parameters:
737: + grid - The grid
738: - alpha - The scalar
740: Level: intermediate
742: .keywords grid, rhs, scale
743: .seealso GridScaleSystemMatrix()
744: @*/
745: int GridScaleRhs(Grid grid, PetscScalar alpha)
746: {
747: int func, op;
751: /* Scale functions */
752: for(func = 0; func < grid->numRhsFuncs; func++) grid->rhsFuncs[func].alpha *= alpha;
753: /* Scale linear operators */
754: for(op = 0; op < grid->numRhsOps; op++) grid->rhsOps[op].alpha *= alpha;
755: PetscLogFlops(grid->numRhsFuncs + grid->numRhsOps);
756: return(0);
757: }
759: #undef __FUNCT__
761: /*@C GridScaleRhsOperator
762: This function multiplies an operator of the Rhs by the a scalar.
764: Collective on Grid
766: Input Parameters:
767: + grid - The grid
768: . alpha - The scalar
769: - index - The unique operator index returned by GridAddRhsOperator()
771: Level: intermediate
773: .keywords grid, rhs, scale, operator
774: .seealso GridScaleSystemMatrix(), GridAddRhsOperator()
775: @*/
776: int GridScaleRhsOperator(Grid grid, PetscScalar alpha, int index)
777: {
780: grid->rhsOps[index].alpha *= alpha;
781: PetscLogFlops(1);
782: return(0);
783: }
785: /*--------------------------------------------- System Matrix Functions ---------------------------------------------*/
786: #undef __FUNCT__
788: /*@
789: GridIsMatrixFree - Signals whether operators are applied matrix-free.
791: Not collective
793: Input Parameter:
794: . grid - The grid
796: Output Parameter:
797: . isMatrixFree - The flag for matrix-free operators
799: Level: intermediate
801: .keywords: grid, operator, matrix-free
802: .seealso: GridSetMatrixFree()
803: @*/
804: int GridIsMatrixFree(Grid grid, PetscTruth *isMatrixFree)
805: {
807: *isMatrixFree = grid->isMatrixFree;
808: return(0);
809: }
811: #undef __FUNCT__
813: /*@
814: GridSetMatrixFree - Determines whether operators are applied matrix-free.
816: Not collective
818: Input Parameters:
819: . grid - The grid
820: . isMatrixFree - The flag for matrix-free operators
822: Level: intermediate
824: .keywords: grid, operator, matrix-free
825: .seealso: GridIsMatrixFree()
826: @*/
827: int GridSetMatrixFree(Grid grid, PetscTruth isMatrixFree)
828: {
830: grid->isMatrixFree = isMatrixFree;
831: return(0);
832: }
834: #undef __FUNCT__
836: /*@
837: GridGetNumRhsOperators - Returns the number of rhs operators for the problem defined on the grid.
839: Not collective
841: Input Parameter:
842: . grid - The grid
844: Output Parameter:
845: . numOps - The number of operators
847: Level: intermediate
849: .keywords: grid, operator
850: .seealso: GridGetRhsOperatorStart()
851: @*/
852: int GridGetNumRhsOperators(Grid grid, int *numOps)
853: {
857: *numOps = grid->numRhsOps;
858: return(0);
859: }
861: #undef __FUNCT__
863: /*@
864: GridGetNumMatOperators - Returns the number of matrix operators for the problem defined on the grid.
866: Not collective
868: Input Parameter:
869: . grid - The grid
871: Output Parameter:
872: . numOps - The number of operators
874: Level: intermediate
876: .keywords: grid, operator
877: .seealso: GridGetMatOperatorStart()
878: @*/
879: int GridGetNumMatOperators(Grid grid, int *numOps)
880: {
884: *numOps = grid->numMatOps;
885: return(0);
886: }
888: #undef __FUNCT__
890: /*@
891: GridGetMatOperatorStart - Retrieves information about the first matrix operator
892: for the problem defined on the grid.
894: Not collective
896: Input Parameter:
897: . grid - The grid
899: Output Parameters:
900: + op - The operator
901: . sField - The shape function field
902: . tField - The test function field
903: . alpha - The multiplier of this operator
904: - isALE - The flag for ALE operators
906: Level: intermediate
908: .keywords: grid, operator
909: .seealso: GridGetMatOperatorNext(), GridGetNumMatOperators()
910: @*/
911: int GridGetMatOperatorStart(Grid grid, int *op, int *sField, int *tField, PetscScalar *alpha, PetscTruth *isALE)
912: {
917: grid->activeMatOp = -1;
918: GridGetMatOperatorNext(grid, op, sField, tField, alpha, isALE);
919: return(0);
920: }
922: #undef __FUNCT__
924: /*@
925: GridGetMatOperatorNext - Retrieves information about the next matrix operator
926: for the problem defined on the grid.
928: Not collective
930: Input Parameter:
931: . grid - The grid
933: Output Parameters:
934: + op - The operator
935: . sField - The shape function field
936: . tField - The test function field
937: . alpha - The multiplier of this operator
938: - isALE - The flag for ALE operators
940: Level: intermediate
942: .keywords: grid, operator
943: .seealso: GridMatGetOperatorStart(), GridGetNumMatOperators()
944: @*/
945: int GridGetMatOperatorNext(Grid grid, int *op, int *sField, int *tField, PetscScalar *alpha, PetscTruth *isALE)
946: {
949: grid->activeMatOp++;
950: if (grid->activeMatOp >= grid->numMatOps) {
951: if (op != PETSC_NULL) *op = -1;
952: if (sField != PETSC_NULL) *sField = -1;
953: if (tField != PETSC_NULL) *tField = -1;
954: if (alpha != PETSC_NULL) *alpha = 0.0;
955: if (isALE != PETSC_NULL) *isALE = PETSC_FALSE;
956: grid->activeMatOp = -1;
957: return(0);
958: }
959: if (op != PETSC_NULL) {
961: *op = grid->matOps[grid->activeMatOp].op;
962: }
963: if (sField != PETSC_NULL) {
965: *sField = grid->matOps[grid->activeMatOp].field;
966: }
967: if (tField != PETSC_NULL) {
969: *tField = grid->fields[grid->matOps[grid->activeMatOp].field].disc->operators[grid->matOps[grid->activeMatOp].op]->test->field;
970: }
971: if (alpha != PETSC_NULL) {
973: *alpha = grid->matOps[grid->activeMatOp].alpha;
974: }
975: if (isALE != PETSC_NULL) {
977: *isALE = grid->matOps[grid->activeMatOp].isALE;
978: }
979: return(0);
980: }
982: #undef __FUNCT__
984: /*@C GridSetMatOperator
985: This function sets the operator to use in the system matrix for the given field.
987: Collective on Grid
989: Input Parameters:
990: + grid - The grid
991: . sField - The field for shape functions
992: . tField - The field for test functions
993: . op - The operator
994: . alpha - A scalar multiple for the operator
995: - isALE - A flag for ALE operators
997: Output Parameter:
998: . index - [Optional] A unique indentifier for the operator
1000: Level: beginner
1002: .keywords active field
1003: .seealso GridAddMatOperator(), GridRemoveMatOperator(), GridAddRhsOperator()
1004: @*/
1005: int GridSetMatOperator(Grid grid, int sField, int tField, int op, PetscScalar alpha, PetscTruth isALE, int *index)
1006: {
1011: grid->numMatOps = 0;
1012: GridAddMatOperator(grid, sField, tField, op, alpha, isALE, index);
1013: return(0);
1014: }
1016: #undef __FUNCT__
1018: /*@C GridAddMatOperator
1019: This function adds the operator to the system matrix for the given field.
1021: Collective on Grid
1023: Input Parameters:
1024: + grid - The grid
1025: . sField - The field for shape functions
1026: . tField - The field for test functions
1027: . op - The operator
1028: . alpha - A scalar multiple for the operator
1029: - isALE - A flag for ALE operators
1031: Output Parameter:
1032: . index - [Optional] A unique indentifier for the operator
1034: Level: beginner
1036: .keywords active field
1037: .seealso GridRemoveMatOperator(), GridSetMatOperator(), GridSetRhsOperator()
1038: @*/
1039: int GridAddMatOperator(Grid grid, int sField, int tField, int op, PetscScalar alpha, PetscTruth isALE, int *index)
1040: {
1041: Discretization sDisc = grid->fields[sField].disc;
1042: Discretization tDisc = grid->fields[tField].disc;
1043: int dim = grid->dim;
1044: GridOp *tempOps;
1045: int ierr;
1049: GridValidField(grid, sField);
1050: GridValidField(grid, tField);
1051: if ((op < 0) || (op >= sDisc->numOps) || (!sDisc->operators[op])) {
1052: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d", op);
1053: }
1054: if ((sDisc->operators[op]->precompInt == PETSC_NULL) &&
1055: (sDisc->operators[op]->opFunc == PETSC_NULL) &&
1056: (sDisc->operators[op]->ALEOpFunc == PETSC_NULL)) {
1057: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d has no associated function", op);
1058: }
1059: if (sDisc->numQuadPoints != tDisc->numQuadPoints) {
1060: SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible quadrature sets");
1061: }
1062: while (grid->numMatOps >= grid->maxMatOps) {
1063: PetscMalloc(grid->maxMatOps*2 * sizeof(GridOp), &tempOps);
1064: PetscLogObjectMemory(grid, grid->maxMatOps * sizeof(GridOp));
1065: PetscMemcpy(tempOps, grid->matOps, grid->maxMatOps * sizeof(GridOp));
1066: PetscFree(grid->matOps);
1067: grid->matOps = tempOps;
1068: grid->maxMatOps *= 2;
1069: }
1070: if (index != PETSC_NULL) {
1072: *index = grid->numMatOps;
1073: }
1074: grid->matOps[grid->numMatOps].op = op;
1075: grid->matOps[grid->numMatOps].field = sField;
1076: grid->matOps[grid->numMatOps].alpha = alpha;
1077: grid->matOps[grid->numMatOps].isALE = isALE;
1078: grid->matOps[grid->numMatOps].nonlinearOp = PETSC_NULL;
1079: grid->numMatOps++;
1080: if (isALE == PETSC_TRUE) grid->ALEActive = PETSC_TRUE;
1081: /* Assign test function field to operator -- Need more sophisticated checking */
1082: switch (op) {
1083: case IDENTITY:
1084: case LAPLACIAN:
1085: if (sDisc->size != tDisc->size) {
1086: SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
1087: }
1088: break;
1089: case GRADIENT:
1090: if (sDisc->comp*dim != tDisc->comp) {
1091: SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
1092: }
1093: break;
1094: case DIVERGENCE:
1095: if (tDisc->comp*dim != sDisc->comp) {
1096: SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
1097: }
1098: break;
1099: }
1100: sDisc->operators[op]->test = tDisc;
1101: return(0);
1102: }
1104: #undef __FUNCT__
1106: /*@C GridRemoveMatOperator
1107: This function removes the operator to the system matrix for the given field.
1109: Collective on Grid
1111: Input Parameters:
1112: + grid - The grid
1113: - index - The operator ID from GridAddMatOperator()
1115: Level: beginner
1117: .keywords active field
1118: .seealso GridAddMatOperator(), GridSetMatOperator(), GridAddRhsOperator()
1119: @*/
1120: int GridRemoveMatOperator(Grid grid, int index) {
1121: int o;
1125: if ((index < 0) || (index >= grid->numMatOps)) {
1126: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid operator ID %d should be in [0, %d)", index, grid->numMatOps);
1127: }
1128: for(o = index; o < grid->numMatOps-1; o++) {
1129: grid->matOps[o] = grid->matOps[o+1];
1130: }
1131: grid->numMatOps--;
1132: return(0);
1133: }
1135: #undef __FUNCT__
1137: /*@C GridRegisterOperator
1138: This function defines a new operator for a certain field
1140: Collective on Grid
1142: Input Parameters:
1143: + grid - The grid
1144: . field - The field on which the operator acts
1145: - func - The function deinfing the operator
1147: Output Parameter:
1148: . newOp - The index of the new operator
1150: Level: intermediate
1152: .keywords grid, operator, register
1153: .seealso GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction
1154: @*/
1155: int GridRegisterOperator(Grid grid, int field, OperatorFunction func, int *newOp)
1156: {
1161: GridValidField(grid, field);
1162: DiscretizationRegisterOperator(grid->fields[field].disc, func, newOp);
1163: return(0);
1164: }
1166: #undef __FUNCT__
1168: /*@C GridRegisterALEOperator
1169: This function defines a new operator for a certain field
1171: Collective on Grid
1173: Input Parameters:
1174: + grid - The grid
1175: . field - The field on which the operator acts
1176: - func - The function deinfing the operator
1178: Output Parameter:
1179: . newOp - The index of the new operator
1181: Level: intermediate
1183: .keywords grid, oeprator, register
1184: .seealso GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction
1185: @*/
1186: int GridRegisterALEOperator(Grid grid, int field, ALEOperatorFunction func, int *newOp)
1187: {
1192: GridValidField(grid, field);
1193: DiscretizationRegisterALEOperator(grid->fields[field].disc, func, newOp);
1194: return(0);
1195: }
1197: #undef __FUNCT__
1199: /*@C GridScaleSystemMatrix
1200: This function multiplies all the current components of the system matrix by the same scalar.
1202: Collective on Grid
1204: Input Parameters:
1205: + grid - The grid
1206: - alpha - The scalar
1208: Level: intermediate
1210: .keywords grid, matrix, scale
1211: .seealso GridScaleRhs()
1212: @*/
1213: int GridScaleSystemMatrix(Grid grid, PetscScalar alpha)
1214: {
1215: int op;
1219: /* Scale linear operators */
1220: for(op = 0; op < grid->numMatOps; op++) grid->matOps[op].alpha *= alpha;
1221: PetscLogFlops(grid->numMatOps);
1222: return(0);
1223: }
1225: #undef __FUNCT__
1227: /*@C GridScaleMatOperator
1228: This function multiplies all the current components of the system matrix by the same scalar.
1230: Collective on Grid
1232: Input Parameters:
1233: + grid - The grid
1234: . alpha - The scalar
1235: - index - The unique operator index returned by GridAddMatOperator()
1237: Level: intermediate
1239: .keywords grid, matrix, scale
1240: .seealso GridScaleRhs(), GridAddMatOperator()
1241: @*/
1242: int GridScaleMatOperator(Grid grid, PetscScalar alpha, int index)
1243: {
1246: grid->matOps[index].alpha *= alpha;
1247: PetscLogFlops(1);
1248: return(0);
1249: }