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, °ree, &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, °ree, &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, °ree, &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, °ree, &support);
1446: return(0);
1447: }