Actual source code: sfcoord.c
1: #include <petsc/private/sfimpl.h>
3: static PetscErrorCode GetBoundingBox_Internal(PetscInt npoints, PetscInt dim, const PetscReal *coords, PetscReal *bbox)
4: {
5: PetscFunctionBegin;
6: for (PetscInt d = 0; d < dim; d++) {
7: bbox[0 * dim + d] = PETSC_MAX_REAL;
8: bbox[1 * dim + d] = PETSC_MIN_REAL;
9: }
10: for (PetscInt i = 0; i < npoints; i++) {
11: for (PetscInt d = 0; d < dim; d++) {
12: bbox[0 * dim + d] = PetscMin(bbox[0 * dim + d], coords[i * dim + d]);
13: bbox[1 * dim + d] = PetscMax(bbox[1 * dim + d], coords[i * dim + d]);
14: }
15: }
16: PetscFunctionReturn(PETSC_SUCCESS);
17: }
19: static PetscBool IntersectBoundingBox_Internal(PetscInt dim, const PetscReal *a, const PetscReal *b, PetscReal tol)
20: {
21: for (PetscInt d = 0; d < dim; d++) {
22: if (a[1 * dim + d] + tol < b[0 * dim + d] || b[1 * dim + d] + tol < a[0 * dim + d]) return PETSC_FALSE;
23: }
24: return PETSC_TRUE;
25: }
27: static PetscBool InBoundingBox_Internal(PetscInt dim, const PetscReal *x, const PetscReal *bbox, PetscReal tol)
28: {
29: for (PetscInt d = 0; d < dim; d++) {
30: if (x[d] + tol < bbox[0 * dim + d] || bbox[1 * dim + d] + tol < x[d]) return PETSC_FALSE;
31: }
32: return PETSC_TRUE;
33: }
35: /*@
36: PetscSFSetGraphFromCoordinates - Create SF by fuzzy matching leaf coordinates to root coordinates
38: Collective
40: Input Parameters:
41: + sf - PetscSF to set graph on
42: . nroots - number of root coordinates
43: . nleaves - number of leaf coordinates
44: . dim - spatial dimension of coordinates
45: . tol - positive tolerance for matching
46: . rootcoords - array of root coordinates in which root i component d is [i*dim+d]
47: - leafcoords - array of root coordinates in which leaf i component d is [i*dim+d]
49: Notes:
50: The tolerance typically represents the rounding error incurred by numerically computing coordinates via
51: possibly-different procedures. Passing anything from `PETSC_SMALL` to `100 * PETSC_MACHINE_EPSILON` is appropriate for
52: most use cases.
54: Example:
55: As a motivating example, consider fluid flow in the x direction with y (distance from a wall). The spanwise direction,
56: z, has periodic boundary conditions and needs some spanwise length to allow turbulent structures to develop. The
57: distribution is stationary with respect to z, so you want to average turbulence variables (like Reynolds stress) over
58: the z direction. It is complicated in a 3D simulation with arbitrary partitioner to uniquely number the nodes or
59: quadrature point coordinates to average these quantities into a 2D plane where they will be visualized, but it's easy
60: to compute the projection of each 3D point into the 2D plane.
62: Suppose a 2D target mesh and 3D source mesh (logically an extrusion of the 2D, though perhaps not created in that way)
63: are distributed independently on a communicator. Each rank passes its 2D target points as root coordinates and the 2D
64: projection of its 3D source points as leaf coordinates. Calling `PetscSFReduceBegin()`/`PetscSFReduceEnd()` on the
65: result will sum data from the 3D sources to the 2D targets.
67: As a concrete example, consider three MPI ranks with targets (roots)
68: .vb
69: Rank 0: (0, 0), (0, 1)
70: Rank 1: (0.1, 0), (0.1, 1)
71: Rank 2: (0.2, 0), (0.2, 1)
72: .ve
73: Note that targets must be uniquely owned. Suppose also that we identify the following leaf coordinates (perhaps via projection from a 3D space).
74: .vb
75: Rank 0: (0, 0), (0.1, 0), (0, 1), (0.1, 1)
76: Rank 1: (0, 0), (0.1, 0), (0.2, 0), (0, 1), (0.1, 1)
77: Rank 2: (0.1, 0), (0.2, 0), (0.1, 1), (0.2, 1)
78: .ve
79: Leaf coordinates may be repeated, both on a rank and between ranks. This example yields the following `PetscSF` capable of reducing from sources to targets.
80: .vb
81: Roots by rank
82: [0] 0: 0.0000e+00 0.0000e+00 0.0000e+00 1.0000e+00
83: [1] 0: 1.0000e-01 0.0000e+00 1.0000e-01 1.0000e+00
84: [2] 0: 2.0000e-01 0.0000e+00 2.0000e-01 1.0000e+00
85: Leaves by rank
86: [0] 0: 0.0000e+00 0.0000e+00 1.0000e-01 0.0000e+00 0.0000e+00
87: [0] 5: 1.0000e+00 1.0000e-01 1.0000e+00
88: [1] 0: 0.0000e+00 0.0000e+00 1.0000e-01 0.0000e+00 2.0000e-01
89: [1] 5: 0.0000e+00 0.0000e+00 1.0000e+00 1.0000e-01 1.0000e+00
90: [1] 10: 2.0000e-01 1.0000e+00
91: [2] 0: 1.0000e-01 0.0000e+00 2.0000e-01 0.0000e+00 1.0000e-01
92: [2] 5: 1.0000e+00 2.0000e-01 1.0000e+00
93: PetscSF Object: 3 MPI processes
94: type: basic
95: [0] Number of roots=2, leaves=4, remote ranks=2
96: [0] 0 <- (0,0)
97: [0] 1 <- (1,0)
98: [0] 2 <- (0,1)
99: [0] 3 <- (1,1)
100: [1] Number of roots=2, leaves=6, remote ranks=3
101: [1] 0 <- (0,0)
102: [1] 1 <- (1,0)
103: [1] 2 <- (2,0)
104: [1] 3 <- (0,1)
105: [1] 4 <- (1,1)
106: [1] 5 <- (2,1)
107: [2] Number of roots=2, leaves=4, remote ranks=2
108: [2] 0 <- (1,0)
109: [2] 1 <- (2,0)
110: [2] 2 <- (1,1)
111: [2] 3 <- (2,1)
112: .ve
114: Level: advanced
116: .seealso: `PetscSFCreate()`, `PetscSFSetGraph()`, `PetscSFCreateByMatchingIndices()`
117: @*/
118: PetscErrorCode PetscSFSetGraphFromCoordinates(PetscSF sf, PetscInt nroots, PetscInt nleaves, PetscInt dim, PetscReal tol, const PetscReal *rootcoords, const PetscReal *leafcoords)
119: {
120: PetscReal bbox[6], *bboxes, *target_coords;
121: PetscMPIInt size, *ranks_needed, num_ranks;
122: PetscInt *root_sizes, *root_starts;
123: PetscSFNode *premote, *lremote;
124: PetscSF psf;
125: MPI_Datatype unit;
126: MPI_Comm comm;
128: PetscFunctionBegin;
129: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
130: PetscCall(GetBoundingBox_Internal(nroots, dim, rootcoords, bbox));
131: PetscCallMPI(MPI_Comm_size(comm, &size));
132: PetscCall(PetscMalloc1(size * 2 * dim, &bboxes));
133: PetscCallMPI(MPI_Allgather(bbox, 2 * dim, MPIU_REAL, bboxes, 2 * dim, MPIU_REAL, comm));
134: PetscCall(GetBoundingBox_Internal(nleaves, dim, leafcoords, bbox));
135: PetscCall(PetscMalloc1(size, &root_sizes));
136: PetscCallMPI(MPI_Allgather(&nroots, 1, MPIU_INT, root_sizes, 1, MPIU_INT, comm));
138: PetscCall(PetscMalloc2(size, &ranks_needed, size + 1, &root_starts));
139: root_starts[0] = 0;
140: num_ranks = 0;
141: for (PetscMPIInt r = 0; r < size; r++) {
142: if (IntersectBoundingBox_Internal(dim, bbox, &bboxes[2 * dim * r], tol)) {
143: ranks_needed[num_ranks++] = r;
144: root_starts[num_ranks] = root_starts[num_ranks - 1] + root_sizes[r];
145: }
146: }
147: PetscCall(PetscFree(root_sizes));
148: PetscCall(PetscMalloc1(root_starts[num_ranks], &premote));
149: for (PetscInt i = 0; i < num_ranks; i++) {
150: for (PetscInt j = root_starts[i]; j < root_starts[i + 1]; j++) {
151: premote[j].rank = ranks_needed[i];
152: premote[j].index = j - root_starts[i];
153: }
154: }
155: PetscCall(PetscSFCreate(comm, &psf));
156: PetscCall(PetscSFSetGraph(psf, nroots, root_starts[num_ranks], NULL, PETSC_USE_POINTER, premote, PETSC_USE_POINTER));
157: PetscCall(PetscMalloc1(root_starts[num_ranks] * dim, &target_coords));
158: PetscCallMPI(MPI_Type_contiguous(dim, MPIU_REAL, &unit));
159: PetscCallMPI(MPI_Type_commit(&unit));
160: PetscCall(PetscSFBcastBegin(psf, unit, rootcoords, target_coords, MPI_REPLACE));
161: PetscCall(PetscSFBcastEnd(psf, unit, rootcoords, target_coords, MPI_REPLACE));
162: PetscCallMPI(MPI_Type_free(&unit));
163: PetscCall(PetscSFDestroy(&psf));
165: // Condense targets to only those that lie within our bounding box
166: PetscInt num_targets = 0;
167: for (PetscInt i = 0; i < root_starts[num_ranks]; i++) {
168: if (InBoundingBox_Internal(dim, &target_coords[i * dim], bbox, tol)) {
169: premote[num_targets] = premote[i];
170: for (PetscInt d = 0; d < dim; d++) target_coords[num_targets * dim + d] = target_coords[i * dim + d];
171: num_targets++;
172: }
173: }
174: PetscCall(PetscFree(bboxes));
175: PetscCall(PetscFree2(ranks_needed, root_starts));
177: PetscCall(PetscMalloc1(nleaves, &lremote));
178: for (PetscInt i = 0; i < nleaves; i++) {
179: for (PetscInt j = 0; j < num_targets; j++) {
180: PetscReal sum = 0;
181: for (PetscInt d = 0; d < dim; d++) sum += PetscSqr(leafcoords[i * dim + d] - target_coords[j * dim + d]);
182: if (sum < tol * tol) {
183: lremote[i] = premote[j];
184: goto matched;
185: }
186: }
187: switch (dim) {
188: case 1:
189: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "No target found for leaf coordinate %g", (double)leafcoords[i * dim + 0]);
190: case 2:
191: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "No target found for leaf coordinate (%g, %g)", (double)leafcoords[i * dim + 0], (double)leafcoords[i * dim + 1]);
192: case 3:
193: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "No target found for leaf coordinate (%g, %g, %g)", (double)leafcoords[i * dim + 0], (double)leafcoords[i * dim + 1], (double)leafcoords[i * dim + 2]);
194: }
195: matched:;
196: }
197: PetscCall(PetscFree(premote));
198: PetscCall(PetscFree(target_coords));
199: PetscCall(PetscSFSetGraph(sf, nroots, nleaves, NULL, PETSC_USE_POINTER, lremote, PETSC_OWN_POINTER));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }