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, msize;
    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:   PetscCall(PetscMPIIntCast(2 * dim, &msize));
    134:   PetscCallMPI(MPI_Allgather(bbox, msize, MPIU_REAL, bboxes, msize, MPIU_REAL, comm));
    135:   PetscCall(GetBoundingBox_Internal(nleaves, dim, leafcoords, bbox));
    136:   PetscCall(PetscMalloc1(size, &root_sizes));
    137:   PetscCallMPI(MPI_Allgather(&nroots, 1, MPIU_INT, root_sizes, 1, MPIU_INT, comm));

    139:   PetscCall(PetscMalloc2(size, &ranks_needed, size + 1, &root_starts));
    140:   root_starts[0] = 0;
    141:   num_ranks      = 0;
    142:   for (PetscMPIInt r = 0; r < size; r++) {
    143:     if (IntersectBoundingBox_Internal(dim, bbox, &bboxes[2 * dim * r], tol)) {
    144:       ranks_needed[num_ranks++] = r;
    145:       root_starts[num_ranks]    = root_starts[num_ranks - 1] + root_sizes[r];
    146:     }
    147:   }
    148:   PetscCall(PetscFree(root_sizes));
    149:   PetscCall(PetscMalloc1(root_starts[num_ranks], &premote));
    150:   for (PetscInt i = 0; i < num_ranks; i++) {
    151:     for (PetscInt j = root_starts[i]; j < root_starts[i + 1]; j++) {
    152:       premote[j].rank  = ranks_needed[i];
    153:       premote[j].index = j - root_starts[i];
    154:     }
    155:   }
    156:   PetscCall(PetscSFCreate(comm, &psf));
    157:   PetscCall(PetscSFSetGraph(psf, nroots, root_starts[num_ranks], NULL, PETSC_USE_POINTER, premote, PETSC_USE_POINTER));
    158:   PetscCall(PetscMalloc1(root_starts[num_ranks] * dim, &target_coords));
    159:   PetscCall(PetscMPIIntCast(dim, &msize));
    160:   PetscCallMPI(MPI_Type_contiguous(msize, MPIU_REAL, &unit));
    161:   PetscCallMPI(MPI_Type_commit(&unit));
    162:   PetscCall(PetscSFBcastBegin(psf, unit, rootcoords, target_coords, MPI_REPLACE));
    163:   PetscCall(PetscSFBcastEnd(psf, unit, rootcoords, target_coords, MPI_REPLACE));
    164:   PetscCallMPI(MPI_Type_free(&unit));
    165:   PetscCall(PetscSFDestroy(&psf));

    167:   // Condense targets to only those that lie within our bounding box
    168:   PetscInt num_targets = 0;
    169:   for (PetscInt i = 0; i < root_starts[num_ranks]; i++) {
    170:     if (InBoundingBox_Internal(dim, &target_coords[i * dim], bbox, tol)) {
    171:       premote[num_targets] = premote[i];
    172:       for (PetscInt d = 0; d < dim; d++) target_coords[num_targets * dim + d] = target_coords[i * dim + d];
    173:       num_targets++;
    174:     }
    175:   }
    176:   PetscCall(PetscFree(bboxes));
    177:   PetscCall(PetscFree2(ranks_needed, root_starts));

    179:   PetscCall(PetscMalloc1(nleaves, &lremote));
    180:   PetscKDTree tree;
    181:   PetscCount *indices;
    182:   PetscReal  *distances;

    184:   PetscCall(PetscKDTreeCreate(num_targets, dim, target_coords, PETSC_USE_POINTER, PETSC_DETERMINE, &tree));
    185:   PetscCall(PetscMalloc2(nleaves, &indices, nleaves, &distances));
    186:   PetscCall(PetscKDTreeQueryPointsNearestNeighbor(tree, nleaves, leafcoords, tol, indices, distances));
    187:   for (PetscInt i = 0; i < nleaves; i++) {
    188:     if (distances[i] < tol) {
    189:       lremote[i] = premote[indices[i]];
    190:     } else {
    191:       switch (dim) {
    192:       case 1:
    193:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "No target found for leaf coordinate %g", (double)leafcoords[i * dim + 0]);
    194:       case 2:
    195:         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]);
    196:       case 3:
    197:         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]);
    198:       }
    199:     }
    200:   }
    201:   PetscCall(PetscFree2(indices, distances));
    202:   PetscCall(PetscKDTreeDestroy(&tree));
    203:   PetscCall(PetscFree(premote));
    204:   PetscCall(PetscFree(target_coords));
    205:   PetscCall(PetscSFSetGraph(sf, nroots, nleaves, NULL, PETSC_USE_POINTER, lremote, PETSC_OWN_POINTER));
    206:   PetscFunctionReturn(PETSC_SUCCESS);
    207: }