Actual source code: elemvec.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: elemvec.c,v 1.13 2000/07/16 23:20:02 knepley Exp $";
  3: #endif

 5:  #include src/grid/gridimpl.h

  7: /* Logging support */
  8: int ELEMENT_VEC_COOKIE;
  9: int ELEMENT_MAT_COOKIE;

 11: /*---------------------------------------------- Element Matrix Functions -------------------------------------------*/
 12: #undef  __FUNCT__
 14: /*@C ElementMatCreate
 15:   This function creates an element matrix for a finite element calculation.

 17:   Collective on MPI_Comm

 19:   Input Parameters:
 20: + comm    - The communicator to associate with the matrix
 21: . rowSize - The number of test functions per element
 22: - colSize - The number of basis functions per element

 24:   Output Parameter:
 25: . mat     - The element matrix

 27:   Level: beginner

 29: .keywords element matrix, finite element
 30: .seealso ElementMatDestroy(), ElementVecCreate()
 31: @*/
 32: int ElementMatCreate(MPI_Comm comm, int rowSize, int colSize, ElementMat *mat)
 33: {
 34:   ElementMat m;
 35:   int        ierr;

 39:   *mat = PETSC_NULL;
 40: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
 41:   GridInitializePackage(PETSC_NULL);
 42: #endif

 44:   PetscHeaderCreate(m, _ElementMat, int, ELEMENT_MAT_COOKIE, 0, "ElementMat", comm, ElementMatDestroy, 0);
 45:   PetscLogObjectCreate(m);
 46:   m->rowSize       = rowSize;
 47:   m->colSize       = colSize;
 48:   m->reduceRowSize = rowSize;
 49:   m->reduceColSize = colSize;
 50:   PetscMalloc(rowSize*colSize            * sizeof(PetscScalar), &m->array);
 51:   PetscMalloc(rowSize                    * sizeof(int),         &m->rowIndices);
 52:   PetscMalloc(colSize                    * sizeof(int),         &m->colIndices);
 53:   PetscMalloc(colSize                    * sizeof(int),         &m->reduceCols);
 54:   PetscMalloc(rowSize*colSize            * sizeof(PetscScalar), &m->tempArray);
 55:   PetscMalloc(PetscMax(rowSize, colSize) * sizeof(int),         &m->tempIndices);
 56:   PetscLogObjectMemory(m, (rowSize + colSize*2 + PetscMax(rowSize, colSize)) * sizeof(int) + (rowSize*colSize*2) * sizeof(PetscScalar));
 57:   *mat = m;

 59:   return(0);
 60: }

 62: #undef  __FUNCT__
 64: /*@C ElementMatDestroy
 65:   This function destroys an element matrix.

 67:   Not collective

 69:   Input Parameter:
 70: . vec - The element matrix

 72:   Level: beginner

 74: .keywords element matrix, finite element
 75: .seealso ElementMatCreate(), ElementVecDestroy()
 76: @*/
 77: int ElementMatDestroy(ElementMat mat)
 78: {

 83:   if (--mat->refct > 0)
 84:     return(0);
 85:   PetscFree(mat->array);
 86:   PetscFree(mat->rowIndices);
 87:   PetscFree(mat->colIndices);
 88:   PetscFree(mat->reduceCols);
 89:   PetscFree(mat->tempArray);
 90:   PetscFree(mat->tempIndices);
 91:   PetscLogObjectDestroy(mat);
 92:   PetscHeaderDestroy(mat);
 93:   return(0);
 94: }

 96: #undef  __FUNCT__
 98: /*@C ElementMatView
 99:   This function destroys an element matrix.

101:   Not collective

103:   Input Parameter:
104: . vec - The element matrix

106:   Level: beginner

108: .keywords element matrix, finite element
109: .seealso ElementMatCreate(), ElementVecView()
110: @*/
111: int ElementMatView(ElementMat mat, PetscViewer viewer)
112: {
113:   int        row, col;
114:   PetscTruth isascii;
115:   int        ierr;

120:   PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &isascii);
121:   if (isascii == PETSC_TRUE) {
122:     if (mat->reduceColSize > 0) {
123:       PetscViewerASCIIPrintf(viewer, "      %3d", mat->colIndices[0]);
124:       for(col = 1; col < mat->reduceColSize; col++)
125:         PetscViewerASCIIPrintf(viewer, "   %3d", mat->colIndices[col]);
126:       PetscViewerASCIIPrintf(viewer, "n");
127:     }
128:     for(row = 0; row < mat->reduceRowSize; row++) {
129:       PetscViewerASCIIPrintf(viewer, "%3d ", mat->rowIndices[row]);
130:       for(col = 0; col < mat->reduceColSize; col++)
131:         PetscViewerASCIIPrintf(viewer, "%5.2g ", PetscRealPart(mat->array[row*mat->reduceColSize+col]));
132:       PetscViewerASCIIPrintf(viewer, "n");
133:     }
134:   }
135:   return(0);
136: }

138: #undef  __FUNCT__
140: /*@C ElementMatDuplicate
141:   This function duplicates an element matrix.

143:   Collective on ElementMat

145:   Input Parameter:
146: . mat    - The original element matrix

148:   Output Parameter:
149: . newmat - The new element matrix

151:   Level: beginner

153: .keywords element matrix, finite element
154: .seealso GridCalcElementMatIndices()
155: @*/
156: int ElementMatDuplicate(ElementMat mat, ElementMat *newmat)
157: {

163:   ElementMatCreate(mat->comm, mat->rowSize, mat->colSize, newmat);
164:   return(0);
165: }

167: #undef  __FUNCT__
169: /*@C ElementMatDuplicateIndices
170:   This function copies the global matrix indices to another element matrix.

172:   Collective on ElementMat

174:   Input Parameter:
175: . mat    - The original element matrix

177:   Output Parameter:
178: . newmat - The element matrix with updated indices

180:   Level: intermediate

182: .keywords element matrix, finite element
183: .seealso GridCalcElementMatIndices()
184: @*/
185: int ElementMatDuplicateIndices(ElementMat mat, ElementMat newmat)
186: {

192:   if (mat == newmat)
193:     return(0);
194:   if ((mat->rowSize != newmat->rowSize) || (mat->colSize != newmat->colSize)) {
195:     SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible matrix sizes");
196:   }
197:   PetscMemcpy(newmat->rowIndices, mat->rowIndices, mat->rowSize * sizeof(int));
198:   PetscMemcpy(newmat->colIndices, mat->colIndices, mat->colSize * sizeof(int));
199:   return(0);
200: }

202: #undef  __FUNCT__
204: /*@C ElementMatZero
205:   This function sets all the entries in the element matrix to zero.

207:   Collective on ElementMat

209:   Output Parameter:
210: . mat - The element matrix

212:   Level: beginner

214: .keywords element matrix, finite element
215: .seealso ElementMatCreate()
216: @*/
217: int ElementMatZero(ElementMat mat)
218: {

223:   PetscMemzero(mat->array, mat->rowSize*mat->colSize * sizeof(PetscScalar));
224:   return(0);
225: }

227: #undef  __FUNCT__
229: /*@C ElementMatSetValues
230:   This function puts all the entries in the element matrix into
231:   the global matrix.

233:   Collective on GMat

235:   Input Parameters:
236: + mat  - The element matrix
237: - addv - The insertion mode
238:  
239:   Output Parameter:
240: . M - The global matrix

242:   Level: beginner

244: .keywords element matrix, finite element
245: .seealso ElementMatCreate
246: @*/
247: int ElementMatSetValues(ElementMat mat, GMat M, InsertMode addv)
248: {

252:   /* Place entries in the global matrix */
255:   MatSetValues(M, mat->reduceRowSize, mat->rowIndices, mat->reduceColSize, mat->colIndices, mat->array, addv);
256: 
257:   mat->reduceRowSize = mat->rowSize;
258:   mat->reduceColSize = mat->colSize;
259:   return(0);
260: }

262: #undef  __FUNCT__
264: /*@C ElementMatSetDiagonalValues
265:   This function puts only the diagonal entries of the element matrix into a global vector.

267:   Collective on GVec

269:   Input Parameters:
270: + mat  - The element matrix
271: - addv - The insertion mode
272:  
273:   Output Parameter:
274: . d    - The global vector

276:   Level: beginner

278: .keywords element matrix, finite element
279: .seealso ElementMatSetValues
280: @*/
281: int ElementMatSetDiagonalValues(ElementMat mat, GVec d, InsertMode addv)
282: {
283:   int *rows = mat->rowIndices;
284:   int *cols = mat->colIndices;
285:   int  row, col;
286:   int  ierr;

289:   /* Place entries in the global matrix */
292:   for(row = 0; row < mat->reduceRowSize; row++) {
293:     for(col = 0; col < mat->reduceColSize; col++) {
294:       if (rows[row] == cols[col]) {
295:         VecSetValue(d, rows[row], mat->array[row*mat->reduceColSize+col], addv);
296:       }
297:     }
298:   }
299:   mat->reduceRowSize = mat->rowSize;
300:   mat->reduceColSize = mat->colSize;
301:   return(0);
302: }

304: #undef  __FUNCT__
306: /*@C ElementMatGetSize
307:   This function returns the size of the element matrix.

309:   Not collective

311:   Input Parameter:
312: . mat     - The element matrix
313:  
314:   Output Parameters:
315: + rowSize - The number of rows
316: - colSize - The number of rows

318:   Level: intermediate

320: .keywords element matrix, finite element
321: .seealso ElementMatGetReduceSize(), ElementMatCreate()
322: @*/
323: int ElementMatGetSize(ElementMat mat, int *rowSize, int *colSize)
324: {
329:   *rowSize = mat->rowSize;
330:   *colSize = mat->colSize;
331:   return(0);
332: }

334: #undef  __FUNCT__
336: /*@C
337:   ElementMatGetReduceSize - This function returns the size of the element matrix
338:   after all reduction have been carried out..

340:   Not collective

342:   Input Parameter:
343: . mat     - The element matrix
344:  
345:   Output Parameters:
346: + rowSize - The number of rows
347: - colSize - The number of rows

349:   Level: intermediate

351: .keywords: element matrix, finite element
352: .seealso: ElementMatGetSize(), ElementMatCreate()
353: @*/
354: int ElementMatGetReduceSize(ElementMat mat, int *rowSize, int *colSize)
355: {
360:   *rowSize = mat->reduceRowSize;
361:   *colSize = mat->reduceColSize;
362:   return(0);
363: }

365: #undef  __FUNCT__
367: /*@C ElementMatGetIndices
368:   This function returns the index mappings for the element matrix.

370:   Collective on ElementMat

372:   Input Parameter:
373: . mat     - The element matrix
374:  
375:   Output Parameters:
376: + rowIdx - The row mapping
377: - colIdx - The column mapping

379:   Level: intermediate

381: .keywords element matrix, finite element
382: .seealso ElementMatCreate
383: @*/
384: int ElementMatGetIndices(ElementMat mat, int **rowIdx, int **colIdx)
385: {
390:   *rowIdx = mat->rowIndices;
391:   *colIdx = mat->colIndices;
392:   return(0);
393: }

395: #undef  __FUNCT__
397: /*@C ElementMatGetArray
398:   This function allows access directly to the element matrix;

400:   Collective on ElementMat

402:   Input Parameters:
403: . mat   - The element matrix
404:  
405:   Output Parameter:
406: . array - The array of values

408:   Level: intermediate

410: .keywords element matrix, finite element
411: .seealso ElementMatCreate
412: @*/
413: int ElementMatGetArray(ElementMat mat, PetscScalar **array)
414: {
418:   *array = mat->array;
419:   return(0);
420: }

422: /*---------------------------------------------- Element Vector Functions -------------------------------------------*/
423: #undef  __FUNCT__
425: /*@C ElementVecCreate
426:   This function creates an element vector for a finite element calculation.

428:   Collective on MPI_Comm

430:   Input Parameters:
431: + comm - The communicator to associate with the vector
432: - size - The number of local variables per element

434:   Output Parameter:
435: . vec - The element vector

437:   Level: beginner

439: .keywords: element vector
440: .seealso: ElementVecDestroy(), ElementMatCreate()
441: @*/
442: int ElementVecCreate(MPI_Comm comm, int size, ElementVec *vec)
443: {
444:   ElementVec v;
445:   int        ierr;

449:   *vec = PETSC_NULL;
450: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
451:   GridInitializePackage(PETSC_NULL);
452: #endif

454:   PetscHeaderCreate(v, _ElementVec, int, ELEMENT_VEC_COOKIE, 0, "ElementVec", comm, ElementVecDestroy, 0);
455:   PetscLogObjectCreate(v);
456:   v->size       = size;
457:   v->reduceSize = size;
458:   PetscMalloc(size * sizeof(PetscScalar), &v->array);
459:   PetscMalloc(size * sizeof(int),         &v->indices);
460:   PetscLogObjectMemory(v, size * (sizeof(int) + sizeof(PetscScalar)));
461:   *vec = v;
462:   return(0);
463: }

465: #undef  __FUNCT__
467: /*@C ElementVecDestroy
468:   This function destroys an element vector.

470:   Not collective

472:   Input Parameter:
473: . vec - The element vector

475:   Level: beginner

477: .keywords: element vector
478: .seealso: ElementVecCreate(), ElementMatDestroy()
479: @*/
480: int ElementVecDestroy(ElementVec vec)
481: {

486:   if (--vec->refct > 0)
487:     return(0);
488:   PetscFree(vec->array);
489:   PetscFree(vec->indices);
490:   PetscLogObjectDestroy(vec);
491:   PetscHeaderDestroy(vec);
492:   return(0);
493: }

495: #undef  __FUNCT__
497: /*@C ElementVecDuplicate
498:   This function duplicates an element vector.

500:   Collective on ElementVec

502:   Input Parameter:
503: . vec    - The original element vector

505:   Output Parameter:
506: . newvec - The new element vector

508:   Level: beginner

510: .keywords element vector, finite element
511: .seealso GridCalcElementVecIndices()
512: @*/
513: int ElementVecDuplicate(ElementVec vec, ElementVec *newvec)
514: {

520:   ElementVecCreate(vec->comm, vec->size, newvec);
521:   return(0);
522: }

524: #undef  __FUNCT__
526: /*@C ElementVecDuplicateIndices
527:   This function copies the global vector indices to another element vector.

529:   Collective on ElementVec

531:   Input Parameter:
532: . vec    - The original element vector

534:   Output Parameter:
535: . newvec - The element vector with updated indices

537:   Level: intermediate

539: .keywords element vector, finite element
540: .seealso GridCalcElementVecIndices()
541: @*/
542: int ElementVecDuplicateIndices(ElementVec vec, ElementVec newvec)
543: {

549:   if (vec == newvec)
550:     return(0);
551:   if (vec->size != newvec->size) SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible vector sizes");
552:   PetscMemcpy(newvec->indices, vec->indices, vec->size * sizeof(int));
553:   return(0);
554: }

556: #undef  __FUNCT__
558: /*@C ElementVecZero
559:   This function sets all the entries in the element vector to zero.

561:   Collective on ElementVec  

563:   Output Parameter:
564: . vec - The element vector

566:   Level: beginner

568: .keywords: element vector
569: .seealso: ElementVecCreate()
570: @*/
571: int ElementVecZero(ElementVec vec)
572: {

577:   PetscMemzero(vec->array, vec->size * sizeof(PetscScalar));
578:   return(0);
579: }

581: #undef  __FUNCT__
583: /*@C ElementVecSetValues
584:   This function puts all the entries in the element vector into
585:   the global vector.

587:   Collective on GVec

589:   Input Parameters:
590: + vec  - The element vector
591: - addv - The insertion mode
592:  
593:   Output Parameter:
594: . M - The global vector

596:   Level: beginner

598: .keywords: element vector
599: .seealso: ElementVecCreate()
600: @*/
601: int ElementVecSetValues(ElementVec vec, GVec v, InsertMode addv)
602: {

606:   /* Place entries in the global vector */
609:   VecSetValues(v, vec->reduceSize, vec->indices, vec->array, addv);
610:   vec->reduceSize = vec->size;
611:   return(0);
612: }

614: #undef  __FUNCT__
616: /*@C ElementVecGetSize
617:   This function returns the size of the element vector.

619:   Not collective

621:   Input Parameter:
622: . vec  - The element vector
623:  
624:   Output Parameter:
625: + size - The number of rows

627:   Level: intermediate

629: .keywords element vector, finite element
630: .seealso ElementVecGetReduceSize(), ElementVecCreate()
631: @*/
632: int ElementVecGetSize(ElementVec vec, int *size)
633: {
637:   *size = vec->size;
638:   return(0);
639: }

641: #undef  __FUNCT__
643: /*@C
644:   ElementVecGetReduceSize - This function returns the size of the element vector
645:   after all reductions have been carried out.

647:   Not collective

649:   Input Parameter:
650: . vec  - The element vector
651:  
652:   Output Parameter:
653: + size - The number of rows

655:   Level: intermediate

657: .keywords: element vector, finite element
658: .seealso: ElementVecGetSize(), ElementVecCreate()
659: @*/
660: int ElementVecGetReduceSize(ElementVec vec, int *size)
661: {
665:   *size = vec->reduceSize;
666:   return(0);
667: }

669: #undef  __FUNCT__
671: /*@C ElementVecGetIndices
672:   This function returns the index mapping for the element vector.

674:   Collective on ElementVec

676:   Input Parameter:
677: . vec - The element vector
678:  
679:   Output Parameter:
680: . idx - The mapping

682:   Level: intermediate

684: .keywords: element vector
685: .seealso: ElementVecCreate()
686: @*/
687: int ElementVecGetIndices(ElementVec vec, int **idx)
688: {
692:   *idx = vec->indices;
693:   return(0);
694: }

696: #undef  __FUNCT__
698: /*@C ElementVecGetArray
699:   This function allows access directly to the element vector;

701:   Collective on ElementVec

703:   Input Parameter:
704: . vec   - The element vector
705:  
706:   Output Parameter:
707: . array - The array of values

709:   Level: intermediate

711: .keywords element vector, finite element
712: .seealso ElementVecRestoreArray(), ElementVecCreate()
713: @*/
714: int ElementVecGetArray(ElementVec vec, PetscScalar **array)
715: {
719:   *array = vec->array;
720:   return(0);
721: }

723: #undef  __FUNCT__
725: /*@C ElementVecRestoreArray
726:   This function relinquishes access to the element vector;

728:   Collective on ElementVec

730:   Input Parameters:
731: + vec   - The element vector
732: - array - The array of values

734:   Level: intermediate

736: .keywords element vector, finite element
737: .seealso ElementVecGetArray(), ElementVecCreate()
738: @*/
739: int ElementVecRestoreArray(ElementVec vec, PetscScalar **array)
740: {
744:   /* Right now we don't check */
745:   return(0);
746: }

748: /*------------------------------------- Element Vector Index Calculation Functions ----------------------------------*/
749: #undef  __FUNCT__
751: /*@C
752:   GridCalcElementVecIndices - This function calculates the global row indices for a particular element.

754:   Not collective

756:   Input Parameters:
757: + grid - The Grid
758: - elem - The element

760:   Output Parameter:
761: . vec - The element vector

763:   Level: advanced

765: .keywords: element vector, index
766: .seealso: GridCalcLocalElementVecIndices(), GridCalcInterpolationElementVecIndices(), GridCalcBoundaryElementVecIndices(),
767:           ElementVecSetValues()
768: @*/
769: int GridCalcElementVecIndices(Grid grid, int elem, ElementVec vec)
770: {

774:   GridCalcGeneralElementVecIndices(grid, elem, grid->order, PETSC_NULL, PETSC_FALSE, vec);
775:   return(0);
776: }

778: #undef  __FUNCT__
780: /*@C
781:   GridCalcLocalElementVecIndices - This function calculates the local row indices for a particular element.

783:   Not collective

785:   Input Parameters:
786: + grid - The Grid
787: - elem - The element

789:   Output Parameter:
790: . vec - The element vector

792:   Level: advanced

794: .keywords: element vector, index, local
795: .seealso: GridCalcElementVecIndices(), GridCalcInterpolationElementVecIndices(), GridCalcBoundaryElementVecIndices(),
796:           ElementVecSetValues()
797: @*/
798: int GridCalcLocalElementVecIndices(Grid grid, int elem, ElementVec vec)
799: {
800:   PetscTruth expConst;
801:   int        ierr;

804:   GridGetExplicitConstraints(grid, &expConst);
805:   if (expConst == PETSC_FALSE) {
806:     GridCalcGeneralElementVecIndices(grid, elem, grid->order,           PETSC_NULL, PETSC_TRUE, vec);
807:   } else {
808:     GridCalcGeneralElementVecIndices(grid, elem, grid->constraintOrder, PETSC_NULL, PETSC_TRUE, vec);
809:   }
810:   return(0);
811: }

813: #undef  __FUNCT__
815: /*@C
816:   GridCalcGeneralElementVecIndices - This function calculates the row and column indices
817:   for a particular element using alternate variable orderings.

819:   Not collective

821:   Input Parameters:
822: + grid           - The Grid
823: . elem           - The element
824: . order          - The global variable ordering
825: . reductionOrder - The ordering on the reduction space, or PETSC_NULL for the default
826: - useLocal       - The flag for local numbering of ghost nodes

828:   Output Parameter:
829: . vec            - The element vector

831:   Level: advanced

833: .keywords: element vector, index, general
834: .seealso: ElementVecSetValues()
835: @*/
836: int GridCalcGeneralElementVecIndices(Grid grid, int elem, VarOrdering order, VarOrdering reductionOrder,
837:                                      PetscTruth useLocal, ElementVec vec)
838: {

845:   if (reductionOrder == PETSC_NULL) {
846:     (*grid->ops->gridcalcelemvecidx)(grid, grid->mesh, elem, order, grid->reduceOrder, useLocal, PETSC_FALSE, vec);
847: 
848:   } else {
850:     (*grid->ops->gridcalcelemvecidx)(grid, grid->mesh, elem, order, reductionOrder, useLocal, PETSC_FALSE, vec);
851: 
852:   }
853:   return(0);
854: }

856: #undef  __FUNCT__
858: /*@C
859:   GridCalcInterpolationElementVecIndices - This function calculates the local row indices for
860:   a particular element using values of the previous mesh and ordering.

862:   Not collective

864:   Input Parameters:
865: + grid  - The Grid
866: . mesh  - The old mesh
867: . elem  - The element
868: - order - The new global variable ordering

870:   Output Parameter:
871: . vec   - The element vector

873:   Level: advanced

875: .keywords: element vector, index, interpolation
876: .seealso: GridCalcElementVecIndices(), GridCalcLocalElementVecIndices(), GridCalcBoundaryElementVecIndices(),
877:           ElementVecSetValues()
878: @*/
879: int GridCalcInterpolationElementVecIndices(Grid grid, Mesh mesh, int elem, VarOrdering order, ElementVec vec)
880: {

888:   (*grid->ops->gridcalcelemvecidx)(grid, mesh, elem, order, grid->reduceOrder, PETSC_TRUE, PETSC_TRUE, vec);
889: 
890:   return(0);
891: }

893: #undef  __FUNCT__
895: /*@C
896:   GridCalcBoundaryElementVecIndices - This function calculates the global row indices for a particular boundary element.

898:   Not collective

900:   Input Parameters:
901: + grid     - The discretization context
902: . bd       - The boundary index
903: . edge     - The boundary element
904: . midNode  - The midnode on an edge, or -1 if none exists
905: . order    - The new global variable ordering
906: - useLocal - The flag for local numbering of ghost nodes

908:   Output Parameter:
909: . vec      - The element vector

911:   Note:
912:   The boundary is specified by its index, not marker. Use
913:   MeshGetBoundaryIndex() to retrieve this from the marker.

915:   Level: advanced

917: .keywords element vector, index, boundary
918: .seealso: GridCalcElementVecIndices(), GridCalcLocalElementVecIndices(), GridCalcInterpolationElementVecIndices(),
919:           ElementVecSetValues()
920: @*/
921: int GridCalcBoundaryElementVecIndices(Grid grid, int bd, int edge, int midNode, VarOrdering order, PetscTruth useLocal,
922:                                       ElementVec vec)
923: {

930:   if ((bd < 0) || (bd >= grid->numBd)) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary index %d", bd);
931:   (*grid->ops->gridcalcboundaryelemvecidx)(grid, bd, edge, midNode, order, grid->reduceOrder, useLocal, vec);
932: 
933:   return(0);
934: }

936: /*----------------------------------------- Element Vector Projection Functions -------------------------------------*/
937: #undef  __FUNCT__
939: /*@C
940:   GridProjectElementVec - This function calculates the global row and column indices
941:   for a particular element using alternate variable orderings, in the space constrained variables.

943:   Not collective

945:   Input Parameters:
946: + grid      - The Grid
947: . mesh      - The mesh
948: . elem      - The element
949: . order     - The global variable ordering
950: - constrain - The flag for constraining (applying P^T), or unconstraining (applying P)

952:   Output Parameter:
953: . vec       - The element vector

955:   Level: advanced

957: .keywords: element vector, constraint, projection
958: .seealso: ElementVecSetValues()
959: @*/
960: int GridProjectElementVec(Grid grid, Mesh mesh, int elem, VarOrdering order, PetscTruth constrain, ElementVec vec)
961: {

967:   (*grid->ops->gridprojectelemvec)(grid, mesh, elem, order, grid->reduceOrder, constrain, PETSC_FALSE, vec);
968: 
969:   return(0);
970: }

972: #undef  __FUNCT__
974: /*@C
975:   GridProjectInterpolationElementVec - This function calculates the global row and column indices
976:   for a particular element using alternate variable orderings, in the space constrained variables,
977:   and uses the previous numbering constructs for interpolation.

979:   Not collective

981:   Input Parameters:
982: + grid      - The Grid
983: . mesh      - The mesh
984: . elem      - The element
985: . order     - The global variable ordering
986: - constrain - The flag for constraining (applying P^T), or unconstraining (applying P)

988:   Output Parameter:
989: . vec       - The element vector

991:   Level: advanced

993: .keywords: element vector, constraint, projection
994: .seealso: ElementVecSetValues()
995: @*/
996: int GridProjectInterpolationElementVec(Grid grid, Mesh mesh, int elem, VarOrdering order, PetscTruth constrain, ElementVec vec)
997: {

1003:   (*grid->ops->gridprojectelemvec)(grid, mesh, elem, order, grid->reduceOrder, constrain, PETSC_TRUE, vec);
1004: 
1005:   return(0);
1006: }

1008: /*------------------------------------- Element Matrix Index Calculation Functions ----------------------------------*/
1009: #undef  __FUNCT__
1011: /*@C
1012:   GridCalcElementMatIndices - This function calculates the global row and column indices for a particular element.

1014:   Not collective

1016:   Input Parameters:
1017: + grid - The Grid
1018: - elem - The element

1020:   Output Parameter:
1021: . mat  - The element matrix

1023:   Level: advanced

1025: .keywords: element matrix, index
1026: .seealso: ElementMatSetValues()
1027: @*/
1028: int GridCalcElementMatIndices(Grid grid, int elem, ElementMat mat)
1029: {

1033:   GridCalcGeneralElementMatIndices(grid, elem, grid->order, grid->order, PETSC_FALSE, mat);
1034:   return(0);
1035: }

1037: #undef  __FUNCT__
1039: /*@C
1040:   GridCalcLocalElementMatIndices - This function calculates the local row and column indices for a particular element.

1042:   Not collective

1044:   Input Parameters:
1045: + grid - The Grid
1046: - elem - The element

1048:   Output Parameter:
1049: . mat  - The element matrix

1051:   Level: advanced

1053: .keywords: element matrix, index, local
1054: .seealso: ElementMatSetValues()
1055: @*/
1056: int GridCalcLocalElementMatIndices(Grid grid, int elem, ElementMat mat)
1057: {

1061:   GridCalcGeneralElementMatIndices(grid, elem, grid->order, grid->order, PETSC_TRUE, mat);
1062:   return(0);
1063: }

1065: #undef  __FUNCT__
1067: /*@C
1068:   GridCalcGeneralElementMatIndices - This function calculates the row and column indices
1069:   for a particular element using alternate variable orderings.

1071:   Not collective

1073:   Input Parameters:
1074: + grid     - The Grid
1075: . elem     - The element
1076: . sOrder   - The global variable ordering for the shape functions 
1077: . tOrder   - The global variable ordering for the test functions 
1078: - useLocal - The flag for local numbering of ghost nodes

1080:   Output Parameter:
1081: . mat      - The element matrix

1083:   Level: advanced

1085: .keywords: element matrix, index, general
1086: .seealso: ElementMatSetValues()
1087: @*/
1088: int GridCalcGeneralElementMatIndices(Grid grid, int elem, VarOrdering sOrder, VarOrdering tOrder, PetscTruth useLocal,
1089:                                      ElementMat mat)
1090: {

1098:   (*grid->ops->gridcalcelemmatidx)(grid, elem, sOrder, tOrder, grid->reduceOrder, useLocal, mat);
1099:   return(0);
1100: }

1102: #undef  __FUNCT__
1104: /*@C
1105:   GridCalcBoundaryElementMatIndices - This function calculates the global row and column indices
1106:   for a particular boundary element using alternate variable orderings.

1108:   Not collective

1110:   Input Parameters:
1111: + grid     - The Grid
1112: . bd       - The boundary index
1113: . edge     - The canonical edge number
1114: . midNode  - [Optional] The canonical node number of the midnode on the edge, or -1 if none exists
1115: . sOrder   - The global variable ordering for the shape functions 
1116: . tOrder   - The global variable ordering for the test functions 
1117: - useLocal - The flag for local numbering of ghost nodes

1119:   Output Parameter:
1120: . mat      - The element matrix

1122:   Note:
1123:   The boundary is specified by its index, not marker. Use
1124:   MeshGetBoundaryIndex() to retrieve this from the marker.

1126:   Level: advanced

1128: .keywords: element matrix, index, boundary
1129: .seealso: ElementMatSetValues()
1130: @*/
1131: int GridCalcBoundaryElementMatIndices(Grid grid, int bd, int edge, int midNode, VarOrdering sOrder, VarOrdering tOrder,
1132:                                       PetscTruth useLocal, ElementMat mat)
1133: {

1141:   if ((bd < 0) || (bd >= grid->numBd)) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary index %d", bd);
1142:   (*grid->ops->gridcalcboundaryelemmatidx)(grid, bd, edge, midNode, sOrder, tOrder, grid->reduceOrder, useLocal, mat);
1143: 
1144:   return(0);
1145: }