Actual source code: meshQuery.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: meshQuery.c,v 1.19 2000/10/17 13:48:57 knepley Exp $";
  3: #endif

  5: /*
  6:      Defines the interface to the mesh query functions
  7: */

 9:  #include src/mesh/meshimpl.h

 11: /*----------------------------------------------- Mesh Query Functions -----------------------------------------------*/
 12: #undef  __FUNCT__
 14: /*@
 15:   MeshSetDimension - This function sets the dimension of the Mesh. This may only be done before a mesh is generated.

 17:   Not collective

 19:   Input Parameters:
 20: + mesh - The mesh
 21: - dim  - The mesh dimension

 23:   Level: intermediate

 25: .keywords: Mesh, dimension
 26: .seealso: MatGetDimension(), PartitionGetDimension()
 27: @*/
 28: int MeshSetDimension(Mesh mesh, int dim)
 29: {
 32:   if ((dim < 0) || (dim > 3)) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Dimension %d is invalid", dim);
 33:   if (mesh->setupcalled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Dimension cannot be set after mesh has been generated");
 34:   mesh->dim = dim;
 35:   return(0);
 36: }

 38: #undef  __FUNCT__
 40: /*@
 41:   MeshGetDimension - This function returns the dimension of the Mesh

 43:   Not collective

 45:   Input Parameter:
 46: . mesh - The mesh

 48:   Output Parameter:
 49: . dim  - The mesh dimension

 51:   Level: intermediate

 53: .keywords: Mesh, dimension
 54: .seealso: MatSetDimension(), PartitionGetDimension()
 55: @*/
 56: int MeshGetDimension(Mesh mesh, int *dim)
 57: {
 61:   *dim = mesh->dim;
 62:   return(0);
 63: }

 65: #undef  __FUNCT__
 67: /*@
 68:   MeshGetInfo - This function returns some information about the mesh.

 70:   Collective on Mesh

 72:   Input Parameter:
 73: . mesh - The mesh

 75:   Output Parameters:
 76: + vertices - The number of vertices
 77: . nodes    - The number of nodes
 78: . edges    - The number of edges
 79: - elements - The number of elements

 81:   Level: intermediate

 83: .keywords: mesh
 84: .seealso: MatGetCorners()
 85: @*/
 86: int MeshGetInfo(Mesh mesh, int *vertices, int *nodes, int *edges, int *elements) {
 89:   if (vertices != PETSC_NULL) *vertices = mesh->numVertices;
 90:   if (nodes    != PETSC_NULL) *nodes    = mesh->numNodes;
 91:   if (edges    != PETSC_NULL) *edges    = mesh->numEdges;
 92:   if (elements != PETSC_NULL) *elements = mesh->numFaces;
 93:   return(0);
 94: }

 96: #undef  __FUNCT__
 98: /*@
 99:   MeshSetNumCorners - This function sets the number of nodes on each element. This may only be done before a mesh is generated.

101:   Not collective

103:   Input Parameters:
104: + mesh       - The mesh
105: - numCorners - The number of nodes on each element

107:   Level: intermediate

109: .keywords: Mesh, dimension
110: .seealso: MatGetNumCorners()
111: @*/
112: int MeshSetNumCorners(Mesh mesh, int numCorners)
113: {
116:   if (numCorners <= 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Number of corners %d is invalid", numCorners);
117:   if (mesh->setupcalled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Corners cannot be set after mesh has been generated");
118:   mesh->numCorners = numCorners;
119:   return(0);
120: }

122: #undef  __FUNCT__
124: /*@
125:   MeshGetNumCorners - This function returns the number of nodes on each element.

127:   Not collective

129:   Input Parameter:
130: . mesh       - The mesh

132:   Output Parameter:
133: . numCorners - The number of nodes on an element

135:   Level: intermediate

137: .keywords: mesh, corner
138: .seealso: MatGetInfo()
139: @*/
140: int MeshGetNumCorners(Mesh mesh, int *numCorners) {
144:   *numCorners = mesh->numCorners;
145:   return(0);
146: }

148: #undef  __FUNCT__
150: /*@
151:   MeshSetBoundingBox - This function sets the bounding box for the mesh. This can only be done before
152:   the mesh is generated.

154:   Not collective

156:   Input Parameters:
157: + mesh                   - The mesh
158: . startX, startY, startZ - The lower-left corner of the box
159: - endX,   endY,   endZ   - The upper-right corner of the box

161:   Level: intermediate

163: .keywords: mesh, bounding box
164: .seealso: MeshGetBoundingBox(), MeshSetLocalBoundingBox(), MeshUpdateBoundingBox()
165: @*/
166: int MeshSetBoundingBox(Mesh mesh, PetscReal startX, PetscReal startY, PetscReal startZ, PetscReal endX, PetscReal endY, PetscReal endZ)
167: {
170:   if ((startX > endX) || (startY > endY) || (startZ > endZ)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "The bounding box is inverted");
171:   if (mesh->setupcalled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Bounding box cannot be set after mesh has been generated");
172:   mesh->startX = startX;
173:   mesh->startY = startY;
174:   mesh->startZ = startZ;
175:   mesh->endX   = endX;
176:   mesh->endY   = endY;
177:   mesh->endZ   = endZ;
178:   mesh->sizeX  = endX - startX;
179:   mesh->sizeY  = endY - startY;
180:   mesh->sizeZ  = endZ - startZ;
181:   return(0);
182: }

184: #undef  __FUNCT__
186: /*@
187:   MeshGetBoundingBox - This function returns the bounding box for the mesh

189:   Not collective

191:   Input Parameter:
192: . mesh - The mesh

194:   Output Parameters:
195: + startX, startY, startZ - The lower-left corner of the box
196: - endX,   endY,   endZ   - The upper-right corner of the box

198:   Level: intermediate

200: .keywords: mesh, bounding box
201: .seealso: MeshSetBoundingBox(), MeshGetLocalBoundingBox(), MeshUpdateBoundingBox()
202: @*/
203: int MeshGetBoundingBox(Mesh mesh, PetscReal *startX, PetscReal *startY, PetscReal *startZ, PetscReal *endX, PetscReal *endY, PetscReal *endZ)
204: {
207:   if (startX != PETSC_NULL) {
209:     *startX = mesh->startX;
210:   }
211:   if (startY != PETSC_NULL) {
213:     *startY = mesh->startY;
214:   }
215:   if (startZ != PETSC_NULL) {
217:     *startZ = mesh->startZ;
218:   }
219:   if (endX   != PETSC_NULL) {
221:     *endX   = mesh->endX;
222:   }
223:   if (endY   != PETSC_NULL) {
225:     *endY   = mesh->endY;
226:   }
227:   if (endZ   != PETSC_NULL) {
229:     *endZ   = mesh->endZ;
230:   }
231:   return(0);
232: }

234: #undef  __FUNCT__
236: /*@
237:   MeshSetLocalBoundingBox - This function sets the local bounding box for the mesh. This can only be done before
238:   the mesh is generated. The local box is the smallest rectangle enclosing the portion of the mesh allocated to
239:   this processor.

241:   Not collective

243:   Input Parameters:
244: + mesh                   - The mesh
245: . startX, startY, startZ - The lower-left corner of the local box
246: - endX,   endY,   endZ   - The upper-right corner of the local box

248:   Level: intermediate

250: .keywords: mesh, bounding box, local
251: .seealso: MeshGetLocalBoundingBox(), MeshSetBoundingBox(), MeshUpdateBoundingBox()
252: @*/
253: int MeshSetLocalBoundingBox(Mesh mesh, PetscReal startX, PetscReal startY, PetscReal startZ, PetscReal endX, PetscReal endY, PetscReal endZ)
254: {
257:   if ((startX > endX) || (startY > endY) || (startZ > endZ)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "The bounding box is inverted");
258:   if (mesh->setupcalled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Bounding box cannot be set after mesh has been generated");
259:   mesh->locStartX = startX;
260:   mesh->locStartY = startY;
261:   mesh->locStartZ = startZ;
262:   mesh->locEndX   = endX;
263:   mesh->locEndY   = endY;
264:   mesh->locEndZ   = endZ;
265:   mesh->locSizeX  = mesh->locEndX - mesh->locStartX;
266:   mesh->locSizeY  = mesh->locEndY - mesh->locStartY;
267:   mesh->locSizeZ  = mesh->locEndZ - mesh->locStartZ;
268:   return(0);
269: }

271: #undef  __FUNCT__
273: /*@
274:   MeshGetLocalBoundingBox - This function returns the local bounding box for the mesh. The local box is the smallest
275:   rectangle enclosing the portion of the mesh allocated to this processor.

277:   Not collective

279:   Input Parameter:
280: . mesh - The mesh

282:   Output Parameters:
283: + startX, startY, startZ - The lower-left corner of the local box
284: - endX,   endY,   endZ   - The upper-right corner of the local box

286:   Level: intermediate

288: .keywords: mesh, bounding box, local
289: .seealso: MeshSetLocalBoundingBox(), MeshGetBoundingBox(), MeshUpdateBoundingBox()
290: @*/
291: int MeshGetLocalBoundingBox(Mesh mesh, PetscReal *startX, PetscReal *startY, PetscReal *startZ, PetscReal *endX, PetscReal *endY, PetscReal *endZ)
292: {
295:   if (startX != PETSC_NULL) {
297:     *startX = mesh->locStartX;
298:   }
299:   if (startY != PETSC_NULL) {
301:     *startY = mesh->locStartY;
302:   }
303:   if (startZ != PETSC_NULL) {
305:     *startZ = mesh->locStartZ;
306:   }
307:   if (endX   != PETSC_NULL) {
309:     *endX   = mesh->locEndX;
310:   }
311:   if (endY   != PETSC_NULL) {
313:     *endY   = mesh->locEndY;
314:   }
315:   if (endZ   != PETSC_NULL) {
317:     *endZ   = mesh->locEndZ;
318:   }
319:   return(0);
320: }

322: #undef  __FUNCT__
324: /*@
325:   MeshUpdateBoundingBox - This function updates the bounding box for the mesh based upon the current geometry.

327:   Not collective

329:   Input Parameter:
330: . mesh - The mesh

332:   Level: intermediate

334: .keywords: mesh, bounding box, update
335: .seealso: MeshSetBoundingBox(), MeshGetBoundingBox()
336: @*/
337: int MeshUpdateBoundingBox(Mesh mesh)
338: {
339:   MPI_Comm   comm;
340:   PetscTruth isPeriodic, isSetup;
341:   PetscReal  locStartX, locStartY, locStartZ;
342:   PetscReal  locEndX,   locEndY,   locEndZ;
343:   PetscReal  startX, startY, startZ;
344:   PetscReal  endX,   endY,   endZ;
345:   int        ierr;

349:   PetscObjectGetComm((PetscObject) mesh, &comm);
350:   MeshIsPeriodic(mesh, &isPeriodic);
351:   if (isPeriodic == PETSC_TRUE) return(0);
352:   /* Calculate the local bounding box */
353:   (*mesh->ops->updateboundingbox)(mesh);
354:   /* Calculate global bounding box */
355:   MeshGetLocalBoundingBox(mesh, &locStartX, &locStartY, &locStartZ, &locEndX, &locEndY, &locEndZ);
356:   MPI_Allreduce(&locStartX, &startX, 1, MPIU_REAL, MPI_MIN, comm);
357:   MPI_Allreduce(&locEndX,   &endX,   1, MPIU_REAL, MPI_MAX, comm);
358:   MPI_Allreduce(&locStartY, &startY, 1, MPIU_REAL, MPI_MIN, comm);
359:   MPI_Allreduce(&locEndY,   &endY,   1, MPIU_REAL, MPI_MAX, comm);
360:   MPI_Allreduce(&locStartZ, &startZ, 1, MPIU_REAL, MPI_MIN, comm);
361:   MPI_Allreduce(&locEndZ,   &endZ,   1, MPIU_REAL, MPI_MAX, comm);
362:   /* Pretend that the mesh is not setup */
363:   isSetup           = mesh->setupcalled;
364:   mesh->setupcalled = PETSC_FALSE;
365:   MeshSetBoundingBox(mesh, startX, startY, startZ, endX, endY, endZ);
366:   mesh->setupcalled = isSetup;
367:   return(0);
368: }

370: #undef  __FUNCT__
372: /*@
373:   MeshGetPartition - Returns the Partition object for this Mesh.

375:   Not collective

377:   Input Parameter:
378: . mesh - The mesh

380:   Output Parameter:
381: . part - The partition

383:   Level: intermediate

385: .keywords: Mesh, Partition, get
386: .seealso: PartitionGetMesh()
387: @*/
388: int MeshGetPartition(Mesh mesh, Partition *part)
389: {
393:   *part = mesh->part;
394:   return(0);
395: }
396: #undef  __FUNCT__
398: /*@
399:   MeshSetMovement - This function sets the dimension of the Mesh. This may only be done before a mesh is generated.

401:   Not collective

403:   Input Parameters:
404: + mesh - The mesh
405: - dim  - The mesh dimension

407:   Level: intermediate

409: .keywords: Mesh, dimension
410: .seealso: MeshGetMovement(), MeshGetMover()
411: @*/
412: int MeshSetMovement(Mesh mesh, PetscTruth isMoving)
413: {
416:   mesh->isMoving = isMoving;
417:   return(0);
418: }

420: #undef  __FUNCT__
422: /*@
423:   MeshGetMovement - This function determines whether the Mesh is moving.

425:   Not collective

427:   Input Parameter:
428: . mesh     - The mesh

430:   Output Parameter:
431: . isMoving - The flag for mesh movement

433:   Level: intermediate

435: .keywords: Mesh, movement
436: .seealso: MeshSetMovement(), MeshGetMover()
437: @*/
438: int MeshGetMovement(Mesh mesh, PetscTruth *isMoving)
439: {
443:   *isMoving = mesh->isMoving;
444:   return(0);
445: }

447: #undef  __FUNCT__
449: /*@
450:   MeshGetMaxDegree - This function returns the maximum degree of any node in the mesh.

452:   Not collective

454:   Input Parameter:
455: . mesh      - The mesh

457:   Output Parameter:
458: . maxDegree - The maximum degree of any node

460:   Level: intermediate

462: .keywords: Mesh, degree
463: .seealso: MeshSetFromOptions()
464: @*/
465: int MeshGetMaxDegree(Mesh mesh, int *maxDegree)
466: {
470:   *maxDegree = mesh->maxDegree;
471:   return(0);
472: }

474: /*------------------------------------------- Mesh Boundary Query Functions ------------------------------------------*/
475: #undef  __FUNCT__
477: /*@
478:   MeshGetNumBoundaries - Returns the number of boundaries in the mesh.

480:   Not collective

482:   Input Parameter:
483: . mesh  - The mesh

485:   Output Parameter:
486: . numBd - The number boundaries

488:   Level: intermediate

490: .keywords: mesh, boundary
491: .seealso: MeshGetBoundaryStart(), MeshGetBoundarySize()
492: @*/
493: int MeshGetNumBoundaries(Mesh mesh, int *numBd) {
497:   *numBd = mesh->numBd;
498:   return(0);
499: }

501: #undef  __FUNCT__
503: /*@
504:   MeshGetBoundarySize - Returns the number of nodes on the specified boundary.

506:   Not collective

508:   Input Parameters:
509: + mesh     - The mesh
510: - boundary - The boundary marker

512:   Output Parameter:
513: . size     - The number of nodes on the boundary

515:   Level: intermediate

517: .keywords: mesh, boundary
518: .seealso: MeshGetBoundaryStart()
519: @*/
520: int MeshGetBoundarySize(Mesh mesh, int boundary, int *size)
521: {

527:   (*mesh->ops->getboundarysize)(mesh, boundary, size);
528:   return(0);
529: }

531: #undef  __FUNCT__
533: /*@
534:   MeshGetBoundaryIndex - Returns the index of the specified boundary.

536:   Not collective

538:   Input Parameters:
539: + mesh     - The mesh
540: - boundary - The boundary marker

542:   Output Parameter:
543: . index - The index of the boundary

545:   Level: intermediate

547: .keywords: mesh, boundary
548: .seealso: MeshGetBoundarySize()
549: @*/
550: int MeshGetBoundaryIndex(Mesh mesh, int boundary, int *index)
551: {

557:   (*mesh->ops->getboundaryindex)(mesh, boundary, index);
558:   return(0);
559: }

561: #undef  __FUNCT__
563: /*@
564:   MeshGetBoundaryStart - Retrieves the canonical node number of the first node
565:   on the specified boundary.

567:   Not collective

569:   Input Parameters:
570: + mesh     - The mesh
571: . boundary - The boundary marker
572: - ghost    - Flag for including ghost nodes

574:   Output Parameter:
575: . node     - The canonical node number

577:   Note:
578: $ This is typically used as an iteration construct with MeshGetBoundaryNext(),
579: $ for example:
580: $
581: $ MeshGetBoundaryStart(mesh, OUTER_BD, &node, PETSC_FALSE);
582: $ while (node >= 0) {
583: $   PetscPrintf(PETSC_COMM_SELF, "I have boundary node %dn", node);
584: $   MeshGetBoundaryNext(mesh, OUTER_BD, &node, PETSC_FALSE);
585: $ }

587:   Level: intermediate

589: .keywords: mesh, boundary, iterator
590: .seealso: MeshGetBoundaryNext()
591: @*/
592: int MeshGetBoundaryStart(Mesh mesh, int boundary, PetscTruth ghost, int *node)
593: {

599:   (*mesh->ops->getboundarystart)(mesh, boundary, ghost, node);
600:   return(0);
601: }

603: #undef  __FUNCT__
605: /*@
606:   MeshGetBoundaryNext - Retrieves the canonical node number of the next node
607:   on the specified boundary, with -1 indicating boundary termination.

609:   Not collective

611:   Input Parameters:
612: + mesh     - The mesh
613: . boundary - The boundary marker
614: - ghost    - Flag for including ghost nodes

616:   Output Parameter:
617: . node     - The canonical node number

619:   Note:
620: $ This is typically used as an iteration construct with MeshGetBoundaryStart(),
621: $ for example:
622: $
623: $ MeshGetBoundaryStart(mesh, OUTER_BD, &node, PETSC_FALSE);
624: $ while (node >= 0) {
625: $   PetscPrintf(PETSC_COMM_SELF, "I have boundary node %dn", node);
626: $   MeshGetBoundaryNext(mesh, OUTER_BD, &node, PETSC_FALSE);
627: $ }

629:   Also, this only returns nodes which lie in the given domain, thus the above
630:   loop would run in parallel on all processors, returning different nodes on each.

632:   Level: intermediate

634: .keywords: mesh, boundary, iterator
635: .seealso: MeshGetBoundaryStart()
636: @*/
637: int MeshGetBoundaryNext(Mesh mesh, int boundary, PetscTruth ghost, int *node)
638: {

644:   (*mesh->ops->getboundarynext)(mesh, boundary, ghost, node);
645:   return(0);
646: }

648: #undef  __FUNCT__
650: /*@
651:   MeshGetActiveBoundary - Retrieves the boundary marker for the boundary
652:   last iterated over, or -1 if no iteration has taken place.

654:   Not collective

656:   Input Parameter:
657: . mesh     - The mesh

659:   Output Parameter:
660: . boundary - The boundary marker

662:   Level: advanced

664: .keywords: grid, boundary, iterator
665: .seealso: GridGetBoundaryNext(), MeshGetBoundaryStart()
666: @*/
667: int MeshGetActiveBoundary(Mesh mesh, int *boundary)
668: {

674:   (*mesh->ops->getactiveboundary)(mesh, boundary);
675:   return(0);
676: }

678: /*--------------------------------------------- Mesh Node Query Functions --------------------------------------------*/
679: #undef  __FUNCT__
681: /*@
682:   MeshGetNodeBoundary - Returns the boundary on which the node lies, or 0 for interior nodes.

684:   Not collective

686:   Input Parameters:
687: + mesh  - The mesh
688: - node  - The node

690:   Output Parameter:
691: . bd    - The boundary

693:   Level: intermediate

695: .keywords: mesh, boundary, node
696: .seealso: MeshGetBoundarySize()
697: @*/
698: int MeshGetNodeBoundary(Mesh mesh, int node, int *bd)
699: {

705:   (*mesh->ops->getnodeboundary)(mesh, node, bd);
706:   return(0);
707: }

709: #undef  __FUNCT__
711: /*@
712:   MeshNodeIsVertex - This function determines whether or not a node is a vertex.

714:   Not collective

716:   Input Parameters:
717: + mesh     - The mesh
718: - node     - The node

720:   Output Parameter:
721: . isVertex - The flag for a vertex

723:   Note:
724:   Vertices are nodes which are connected to edges.

726:   Level: intermediate

728: .keywords: mesh, node, vertex
729: .seealso: MeshGetNodeCoords()
730: @*/
731: int MeshNodeIsVertex(Mesh mesh, int node, PetscTruth *isVertex)
732: {

738:   (*mesh->ops->nodeisvertex)(mesh, node, isVertex);
739:   return(0);
740: }

742: #undef  __FUNCT__
744: /*@
745:   MeshGetNodeCoords - Retrieves the coordinates of a selected node

747:   Input Parameters:
748: + mesh  - The mesh
749: - node  - The local node number

751:   Output Parameters:
752: . x,y,z - The coordinates

754:   Level: intermediate

756: .keywords: mesh, hole, coordinates
757: .seealso: MeshGetBoundaryStart()
758: @*/
759: int MeshGetNodeCoords(Mesh mesh, int node, double *x, double *y, double *z)
760: {

768:   (*mesh->ops->getnodecoords)(mesh, node, x, y, z);
769:   return(0);
770: }

772: #undef  __FUNCT__
774: /*@
775:   MeshSetNodeCoords - Sets the coordinates of a selected node

777:   Collective on Mesh

779:   Input Parameters:
780: + mesh  - The mesh
781: . node  - The local node number
782: - x,y,z - The coordinates

784:   Level: intermediate

786: .keywords: mesh, node, coordinates
787: .seealso: MeshGetBoundaryStart()
788: @*/
789: int MeshSetNodeCoords(Mesh mesh, int node, double x, double y, double z)
790: {

795:   (*mesh->ops->setnodecoords)(mesh, node, x, y, z);
796:   return(0);
797: }

799: #undef  __FUNCT__
801: /*@
802:   MeshGetNodeCoordsSaved - Retrieves the coordinates of a selected node
803:   from those saved with MeshSaveMesh().

805:   Input Parameters:
806: + mesh  - The mesh
807: - node  - The canonical node number

809:   Output Parameters:
810: . x,y,z - The coordinates

812:   Level: intermediate

814: .keywords: mesh, node, coordinates
815: .seealso: MeshGetNodeCoords(), MeshSaveMesh()
816: @*/
817: int MeshGetNodeCoordsSaved(Mesh mesh, int node, double *x, double *y, double *z)
818: {

826:   (*mesh->ops->getnodecoordssaved)(mesh, node, x, y, z);
827:   return(0);
828: }

830: #undef  __FUNCT__
832: /*@
833:   MeshGetNearestNode - This function returns the node nearest to
834:   the given point, or -1 if the closest node is not contained in
835:   the local mesh.

837:   Not collective

839:   Input Parameters:
840: + mesh    - The mesh
841: . x,y,z   - The node coordinates
842: - outside - A flag to allow points outside the domain

844:   Output Parameter:
845: . node    - The nearest node

847:   Note:
848:   The outside flag allows points outside the domain to be tested. If this flag
849:   is PETSC_FALSE, and (x,y,z) does not lie in the global domain, then an error
850:   will result.

852:   Level: beginner

854: .keywords: mesh, node, point location
855: .seealso MeshLocatePoint()
856: @*/
857: int MeshGetNearestNode(Mesh mesh, double x, double y, double z, PetscTruth outside, int *node)
858: {

864:   (*mesh->ops->nearestnode)(mesh, x, y, z, outside, node);
865:   return(0);
866: }

868: #undef  __FUNCT__
870: /*@
871:   MeshGetNearestBdNode - This function returns the boundary node nearest to
872:   the given point, or -1 if the closest node is not contained in the local mesh.

874:   Not collective

876:   Input Parameters:
877: + mesh  - The mesh
878: - x,y,z - The node coordinates

880:   Output Parameter:
881: . node  - The nearest boundary node

883:   Level: beginner

885: .keywords: mesh, node, point location
886: .seealso MeshGetNearestNode(), MeshLocatePoint()
887: @*/
888: int MeshGetNearestBdNode(Mesh mesh, double x, double y, double z, int *node)
889: {
890:   Partition p;
891:   double    minDist = PETSC_MAX;
892:   double    dist, nx, ny, nz;
893:   int       minNode, globalMinNode, allMinNode;
894:   int       numCorners, startNode, endNode;
895:   int       elem, corner, nearNode, bd;
896:   int       ierr;

901:   MeshGetNumCorners(mesh, &numCorners);
902:   MeshGetPartition(mesh, &p);
903:   PartitionGetStartNode(p, &startNode);
904:   PartitionGetEndNode(p, &endNode);
905:   MeshLocatePoint(mesh, x, y, z, &elem);
906:   if (elem >= 0) {
907:     /* Find first boundary node */
908:     for(corner = 0; corner < numCorners; corner++) {
909:       MeshGetNodeFromElement(mesh, elem, corner, &minNode);
910:       MeshGetNodeBoundary(mesh, minNode, &bd);
911:       if (bd != 0) {
912:         MeshGetNodeCoords(mesh, minNode, &nx, &ny, &nz);
913:         minDist = PetscSqr(MeshPeriodicDiffX(mesh, nx - x)) + PetscSqr(MeshPeriodicDiffY(mesh, ny - y));
914:         break;
915:       }
916:     }
917:     if (corner == numCorners) SETERRQ1(PETSC_ERR_ARG_WRONG, "Element %d has no node on a boundary", elem);
918:     /* Find closest node with field defined */
919:     for(corner = 1; corner < numCorners; corner++) {
920:       MeshGetNodeFromElement(mesh, elem, corner, &nearNode);
921:       MeshGetNodeCoords(mesh, nearNode, &nx, &ny, &nz);
922:       MeshGetNodeBoundary(mesh, nearNode, &bd);
923:       dist = PetscSqr(nx - x) + PetscSqr(ny - y);
924:       if ((bd != 0) && (dist < minDist)) {
925:         minDist = dist;
926:         minNode = nearNode;
927:       }
928:     }
929:     PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode);
930:   } else {
931:     minNode       = -1;
932:     globalMinNode = -1;
933:   }

935:   /* The minimum node might still be a ghost node */
936:   MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, mesh->comm);
937:   if ((allMinNode >= startNode) && (allMinNode < endNode)) {
938:     *node = allMinNode - startNode;
939:   } else {
940:     *node = -1;
941:   }
942:   if (allMinNode < 0) PetscFunctionReturn(1);
943:   return(0);
944: }

946: #undef  __FUNCT__
948: /*@
949:   MeshGetNodeSupport - This function get the canonical element numbers of all elements
950:   within the support of a given basis function. A call to MeshRestoreNodeSupport() must
951:   be made before another call to this function.

953:   Not collective

955:   Input Parameters:
956: + mesh    - The mesh
957: . node    - The node containing the basis function
958: - elem    - An element containing the node, -1 for a search

960:   Output Parameters:
961: + degree  - The degree of the node
962: - support - A list of the elements in the support

964:   Note:
965:   This function currently only returns the elements containing
966:   any given node, so some basis functions will have a wider
967:   support than this definition.

969:   Level: intermediate

971: .keywords: mesh, support
972: .seealso: MeshRestoreNodeSupport()
973: @*/
974: int MeshGetNodeSupport(Mesh mesh, int node, int elem, int *degree, int **support)
975: {

982:   if (mesh->supportTaken == PETSC_TRUE) {
983:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Missing call to MeshRestoreNodeSupport()");
984:   }
985:   mesh->supportTaken = PETSC_TRUE;
986:   (*mesh->ops->getnodesupport)(mesh, node, elem, degree, support);
987:   return(0);
988: }

990: #undef  __FUNCT__
992: /*@
993:   MeshRestoreNodeSupport - This function returns the support array from MeshGetNodeSupport().

995:   Not collective

997:   Input Parameters:
998: + mesh    - The mesh
999: . node    - The node containing the basis function
1000: . elem    - An element containing the node, -1 for a search
1001: . degree  - The degree pf the node
1002: - support - A list of the elements in the support

1004:   Level: intermediate

1006: .keywords: mesh, support
1007: .seealso: MeshGetNodeSupport()
1008: @*/
1009: int MeshRestoreNodeSupport(Mesh mesh, int node, int elem, int *degree, int **support)
1010: {
1013:   if (*support != mesh->support) {
1014:     SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid support argument");
1015:   }
1016:   mesh->supportTaken = PETSC_FALSE;
1017:   return(0);
1018: }

1020: #undef  __FUNCT__
1022: /*@
1023:   MeshGetNodeOrdering - This function gets the AO which was used to reorder the mesh nodes
1024:   before partitioning. This is sometimes used to bandwdith reduce the mesh graph.

1026:   Not collective

1028:   Input Parameter:
1029: . mesh  - The mesh

1031:   Output Parameter:
1032: . order - The node ordering, or PETSC_NULL if none exists

1034:   Level: intermediate

1036: .keywords: Mesh, node, ordering, AO
1037: .seealso: PartitionGetNodeOrdering()
1038: @*/
1039: int MeshGetNodeOrdering(Mesh mesh, AO *order)
1040: {
1044:   *order = mesh->nodeOrdering;
1045:   return(0);
1046: }

1048: /*------------------------------------------- Mesh Element Query Functions -------------------------------------------*/
1049: #undef  __FUNCT__
1051: /*@
1052:   MeshGetElementNeighbor - This function retrieves the local element number of the element, or neighbor, opposing a
1053:   given corner.

1055:   Not collective

1057:   Input Parameters:
1058: + mesh     - The mesh
1059: . elem     - The element
1060: - corner   - The corner, or element node number

1062:   Output Parameter:
1063: . neighbor - The neighbor opposing corner

1065:   Level: intermediate

1067: .keywords: mesh, element, neighbor
1068: .seealso: MeshGetNodeCoords()
1069: @*/
1070: int MeshGetElementNeighbor(Mesh mesh, int elem, int corner, int *neighbor)
1071: {

1077:   (*mesh->ops->getelemneighbor)(mesh, elem, corner, neighbor);
1078:   return(0);
1079: }

1081: #undef  __FUNCT__
1083: /*@
1084:   MeshLocatePoint - This function returns the element containing
1085:   the given point, or -1 if it is not contained in the local mesh.

1087:   Not collective

1089:   Input Parameters:
1090: + mesh           - The mesh
1091: - x,y,z          - The node coordinates

1093:   Output Parameter:
1094: . containingElem - An element containing the node, -1 for failure

1096:   Level: beginner

1098: .keywords: mesh, point location
1099: .seealso MeshGetNearestNode(), MeshGetNodeSupport()
1100: @*/
1101: int MeshLocatePoint(Mesh mesh, double x, double y, double z, int *containingElem)
1102: {

1108:   PetscLogEventBegin(MESH_LocatePoint, mesh, 0, 0, 0);
1109:   (*mesh->ops->locatepoint)(mesh, x, y, z, containingElem);
1110:   PetscLogEventEnd(MESH_LocatePoint, mesh, 0, 0, 0);
1111:   return(0);
1112: }

1114: /*--------------------------------------------- Mesh Hole Query Functions --------------------------------------------*/
1115: #undef  __FUNCT__
1117: /*@
1118:   MeshSetHoleCoords - Sets the coordinates of holes in the mesh

1120:   Collective on Mesh

1122:   Input Parameter:
1123: + mesh   - The mesh
1124: . num    - The number of holes
1125: - coords - The coordinates

1127:   Level: advanced

1129: .keywords: mesh, hole, coordinates
1130: .seealso: MeshGetBoundaryStart()
1131: @*/
1132: int MeshSetHoleCoords(Mesh mesh, int num, Vec coords)
1133: {
1134:   int          dim = mesh->dim;
1135:   PetscScalar *array;
1136:   int          size, vecDim;
1137:   int          h, d;
1138:   int          ierr;

1143:   if (num != mesh->numHoles) {
1144:     mesh->numHoles = num;
1145:     if (mesh->holes != PETSC_NULL) {
1146:       PetscFree(mesh->holes);
1147:     }
1148:     if (mesh->numHoles > 0) {
1149:       PetscMalloc(mesh->numHoles*dim * sizeof(double), &mesh->holes);
1150:     } else {
1151:       mesh->holes = PETSC_NULL;
1152:     }
1153:   }
1154:   VecGetArray(coords, &array);
1155:   VecGetLocalSize(coords, &size);
1156:   if (size == mesh->numHoles*dim) {
1157:     PetscMemcpy(mesh->holes, array, mesh->numHoles*dim * sizeof(double));
1158:   } else if ((size > mesh->numHoles*dim) && (size%mesh->numHoles == 0)) {
1159:     vecDim = (int) (size/mesh->numHoles);
1160:     for(h = 0; h < mesh->numHoles; h++) {
1161:       for(d = 0; d < dim; d++) {
1162:         mesh->holes[h*dim+d] = PetscRealPart(array[h*vecDim+d]);
1163:       }
1164:     }
1165:   } else {
1166:     SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid vector of hole coordinates");
1167:   }
1168:   VecRestoreArray(coords, &array);
1169:   return(0);
1170: }

1172: /*------------------------------------------ Mesh Embedding Query Functions ------------------------------------------*/
1173: #undef  __FUNCT__
1175: /*@
1176:   MeshGetElementFromNode - This function retrieves a local element number from the local
1177:   node number. Notice that this could be nondeterministic.

1179:   Not collective

1181:   Input Parameters:
1182: + mesh - The mesh
1183: - node - The number

1185:   Output Parameter:
1186: . elem - The element

1188:   Level: intermediate

1190: .keywords: mesh, node, element
1191: .seealso: MeshGetNodeFromElement()
1192: @*/
1193: int MeshGetElementFromNode(Mesh mesh, int node, int *elem)
1194: {
1195:   Partition part;
1196:   int       numElements, numCorners;
1197:   int       e, c, n;
1198:   int       ierr;

1203:   MeshGetPartition(mesh, &part);
1204:   PartitionGetNumElements(part, &numElements);
1205:   MeshGetNumCorners(mesh, &numCorners);
1206:   for(e = 0; e < numElements; e++) {
1207:     for(c = 0; c < numCorners; c++){
1208:       MeshGetNodeFromElement(mesh, e, c, &n);
1209:       if (n == node) {
1210:         *elem = e;
1211:         return(0);
1212:       }
1213:     }
1214:   }
1215:   PetscFunctionReturn(-1);
1216: }

1218: #undef  __FUNCT__
1220: /*@
1221:   MeshGetBdElementFromEdge - This function retrieves the local element number of a
1222:   boundary element from the local edge number.

1224:   Not collective

1226:   Input Parameters:
1227: + mesh - The mesh
1228: - edge - The edge

1230:   Output Parameter:
1231: . elem - The element along the boundary

1233:   Level: intermediate

1235: .keywords: mesh, element, edge
1236: .seealso: MeshGetNodeFromElement(), MeshGetNodeCoords()
1237: @*/
1238: int MeshGetBdElementFromEdge(Mesh mesh, int edge, int *elem)
1239: {
1240: #if 0
1241:   int  neighbor;
1242: #endif
1243:   int  degree;
1244:   int *support;
1245:   int  startNode,   endNode;
1246:   int  startCorner, endCorner, numCorners;
1247:   int  sElem, corner, node;
1248:   int  ierr;

1252:   /* This should be a default function, which is only used if the implementation
1253:      does not have an optimized implementation
1254:   */
1255:   /* Search for element containing edge */
1256:   MeshGetNumCorners(mesh, &numCorners);
1257:   MeshGetNodeFromEdge(mesh, edge, 0, &startNode);
1258:   MeshGetNodeFromEdge(mesh, edge, 1, &endNode);
1259:   /* Loop over nodes on each element in the support of the node */
1260:   MeshGetNodeSupport(mesh, startNode, 0, &degree, &support);
1261:   for(sElem = 0; sElem < degree; sElem++) {
1262:     for(corner = 0, startCorner = -1, endCorner = -1; corner < numCorners; corner++) {
1263:       MeshGetNodeFromElement(mesh, support[sElem], corner, &node);
1264:       if (node == startNode) {
1265:         startCorner = corner;
1266:       } else if (node == endNode) {
1267:         endCorner   = corner;
1268:       }
1269:     }
1270:     if ((startCorner >= 0) && (endCorner >= 0)) break;
1271:   }
1272:   if (sElem == degree) SETERRQ1(PETSC_ERR_PLIB, "Edge %d not found in element list", edge);
1273:   *elem = support[sElem];
1274:   MeshRestoreNodeSupport(mesh, startNode, 0, &degree, &support);
1275: #if 0
1276:   MeshGetNeighbor(mesh, *elem, ((startCorner+endCorner)*2)%3, &neighbor);
1277:   if (neighbor != -1)
1278:     SETERRQ2(PETSC_ERR_ARG_WRONG, "Edge is not on a boundary since elem %d has neighbor %d", *elem, neighbor);
1279: #endif
1280:   return(0);
1281: }

1283: #undef  __FUNCT__
1285: /*@
1286:   MeshGetNodeFromElement - This function retrieves the local node number from the local
1287:   element number and the corner.

1289:   Not collective

1291:   Input Parameters:
1292: + mesh   - The mesh
1293: . elem   - The element
1294: - corner - The corner, or element node number

1296:   Output Parameter:
1297: . node   - The local node number

1299:   Level: intermediate

1301: .keywords: mesh, node
1302: .seealso: MeshGetNodeCoords()
1303: @*/
1304: int MeshGetNodeFromElement(Mesh mesh, int elem, int corner, int *node)
1305: {
1306:   int numCorners         = mesh->numCorners;
1307:   int numOverlapElements = mesh->numFaces;

1313:   if (mesh->part != PETSC_NULL) {
1314:     ierr  = PartitionGetNumOverlapElements(mesh->part, &numOverlapElements);
1315:   }
1316:   if ((elem < 0) || (elem >= numOverlapElements)) {
1317:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node %d should be in [0,%d)", elem, numOverlapElements);
1318:   }
1319:   if ((corner < 0) || (corner >= numCorners)) {
1320:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid corner %d should be in [0,%d)", elem, numCorners);
1321:   }
1322:   (*mesh->ops->getnodefromelement)(mesh, elem, corner, node);
1323:   return(0);
1324: }

1326: #undef  __FUNCT__
1328: /*@
1329:   MeshGetNodeFromEdge - This function retrieves the local node number from the local
1330:   edge number and the corner.

1332:   Not collective

1334:   Input Parameters:
1335: + mesh   - The mesh
1336: . edge   - The edge
1337: - corner - The corner, or edge node number (0 or 1)

1339:   Output Parameter:
1340: . node   - The local node number

1342:   Level: intermediate

1344: .keywords: mesh, node, edge
1345: .seealso: MeshGetNodeCoords()
1346: @*/
1347: int MeshGetNodeFromEdge(Mesh mesh, int edge, int corner, int *node)
1348: {
1349:   int numOverlapEdges = mesh->numEdges;

1355:   if ((corner < 0) || (corner >= 2)) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid corner %d should be in [0,2)", edge);
1356:   if (mesh->part != PETSC_NULL) {
1357:     ierr  = PartitionGetNumOverlapEdges(mesh->part, &numOverlapEdges);
1358:   }
1359:   if ((edge < 0) || (edge >= numOverlapEdges)) {
1360:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid edge %d should be in [0,%d)", edge, numOverlapEdges);
1361:   }
1362:   (*mesh->ops->getnodefromedge)(mesh, edge, corner, node);
1363:   return(0);
1364: }

1366: #undef  __FUNCT__
1368: /*@
1369:   MeshGetMidnodeFromEdge - This function retrieves the local node number of the midnode from the local edge number.

1371:   Not collective

1373:   Input Parameters:
1374: + mesh    - The mesh
1375: - edge    - The edge

1377:   Output Parameter:
1378: . midnode - The local node number of the midnode

1380:   Level: intermediate

1382: .keywords: mesh, midnode, edge
1383: .seealso: MeshGetNodeFromElement(), MeshGetNodeCoords()
1384: @*/
1385: int MeshGetMidnodeFromEdge(Mesh mesh, int edge, int *midnode)
1386: {
1387:   int  degree;
1388:   int *support;
1389:   int  startNode,   endNode;
1390:   int  numCorners;
1391:   int  startCorner = -1;
1392:   int  endCorner   = -1;
1393:   int  elem, corner, node;
1394:   int  ierr;

1398:   /* This should be a default function, which is only used if the implementation
1399:      does not have an optimized implementation
1400:   */
1401:   MeshGetNumCorners(mesh, &numCorners);
1402:   if (numCorners == 3) return(0);
1403:   if (numCorners != 6) SETERRQ1(PETSC_ERR_ARG_WRONG, "Cannot find midnode for %d corners", numCorners);

1405:   /* Search for midnode on edge */
1406:   MeshGetNodeFromEdge(mesh, edge, 0, &startNode);
1407:   MeshGetNodeFromEdge(mesh, edge, 1, &endNode);
1408:   /* Loop over nodes on each element in the support of the node */
1409:   MeshGetNodeSupport(mesh, startNode, 0, &degree, &support);
1410:   for(elem = 0; elem < degree; elem++) {
1411:     for(corner = 0, startCorner = -1, endCorner = -1; corner < numCorners; corner++) {
1412:       MeshGetNodeFromElement(mesh, support[elem], corner, &node);
1413:       if (node == startNode) {
1414:         startCorner = corner;
1415:       } else if (node == endNode) {
1416:         endCorner   = corner;
1417:       }
1418:     }
1419:     if ((startCorner >= 0) && (endCorner >= 0)) break;
1420:   }
1421:   if (elem == degree) {
1422:     if (startCorner < 0) SETERRQ1(PETSC_ERR_PLIB, "Invalid support of node %d did not contain node", startNode);
1423:     if (endCorner   < 0) {
1424:       for(elem = 0; elem < mesh->numFaces; elem++) {
1425:         for(corner = 0, startCorner = -1, endCorner = -1; corner < numCorners; corner++) {
1426:           MeshGetNodeFromElement(mesh, elem, corner, &node);
1427:           if (node == startNode) {
1428:             startCorner = corner;
1429:           } else if (node == endNode) {
1430:             endCorner   = corner;
1431:           }
1432:         }
1433:         if ((startCorner >= 0) && (endCorner >= 0)) break;
1434:       }
1435:       if (elem == mesh->numFaces) {
1436:         SETERRQ3(PETSC_ERR_PLIB, "Edge %d (%d,%d) not found in element list", edge, startNode, endNode);
1437:       } else {
1438:         SETERRQ5(PETSC_ERR_PLIB, "Edge %d (%d,%d) found in element %d not in support list of node %d",
1439:                  edge, startNode, endNode, elem, startNode);
1440:       }
1441:     }
1442:   }

1444:   MeshGetNodeFromElement(mesh, support[elem], ((startCorner + endCorner)*2)%3 + 3, midnode);
1445:   MeshRestoreNodeSupport(mesh, startNode, 0, &degree, &support);
1446:   return(0);
1447: }