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