Actual source code: ex17.c

  1: static char help[] = "Tests for point location\n\n";

  3: #include <petscsf.h>
  4: #include <petscdmplex.h>

  6: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
  7: {
  9:   DMCreate(comm, dm);
 10:   DMSetType(*dm, DMPLEX);
 11:   DMSetFromOptions(*dm);
 12:   DMViewFromOptions(*dm, NULL, "-dm_view");
 13:   return 0;
 14: }

 16: static PetscErrorCode TestLocation(DM dm)
 17: {
 18:   Vec                points;
 19:   PetscSF            cellSF = NULL;
 20:   const PetscSFNode *cells;
 21:   PetscScalar       *a;
 22:   PetscInt           cdim, n;
 23:   PetscInt           cStart, cEnd, c;

 26:   DMGetCoordinateDim(dm, &cdim);
 27:   DMGetCoordinatesLocalSetUp(dm);
 28:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
 29:   /* Locate all centroids */
 30:   VecCreateSeq(PETSC_COMM_SELF, (cEnd - cStart) * cdim, &points);
 31:   VecSetBlockSize(points, cdim);
 32:   VecGetArray(points, &a);
 33:   for (c = cStart; c < cEnd; ++c) {
 34:     PetscReal centroid[3];
 35:     PetscInt  off = (c - cStart) * cdim, d;

 37:     DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
 38:     for (d = 0; d < cdim; ++d) a[off + d] = centroid[d];
 39:   }
 40:   VecRestoreArray(points, &a);
 41:   DMLocatePoints(dm, points, DM_POINTLOCATION_NONE, &cellSF);
 42:   VecDestroy(&points);
 43:   PetscSFGetGraph(cellSF, NULL, &n, NULL, &cells);
 44:   if (n != (cEnd - cStart)) {
 45:     for (c = 0; c < n; ++c) {
 46:       if (cells[c].index != c + cStart) PetscPrintf(PETSC_COMM_SELF, "Could not locate centroid of cell %" PetscInt_FMT ", error %" PetscInt_FMT "\n", c + cStart, cells[c].index);
 47:     }
 48:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Located %" PetscInt_FMT " points instead of %" PetscInt_FMT, n, cEnd - cStart);
 49:   }
 51:   PetscSFDestroy(&cellSF);
 52:   return 0;
 53: }

 55: int main(int argc, char **argv)
 56: {
 57:   DM dm;

 60:   PetscInitialize(&argc, &argv, NULL, help);
 61:   CreateMesh(PETSC_COMM_WORLD, &dm);
 62:   TestLocation(dm);
 63:   DMDestroy(&dm);
 64:   PetscFinalize();
 65:   return 0;
 66: }

 68: /*TEST

 70:   testset:
 71:     args: -dm_plex_dim 1 -dm_plex_box_faces 10

 73:     test:
 74:       suffix: seg

 76:     test:
 77:       suffix: seg_hash
 78:       args: -dm_refine 2 -dm_plex_hash_location

 80:   testset:
 81:     args: -dm_plex_box_faces 5,5

 83:     test:
 84:       suffix: tri
 85:       requires: triangle

 87:     test:
 88:       suffix: tri_hash
 89:       requires: triangle
 90:       args: -dm_refine 2 -dm_plex_hash_location

 92:     test:
 93:       suffix: quad
 94:       args: -dm_plex_simplex 0

 96:     test:
 97:       suffix: quad_hash
 98:       args: -dm_plex_simplex 0 -dm_refine 2 -dm_plex_hash_location

100:   testset:
101:     args: -dm_plex_dim 3 -dm_plex_box_faces 3,3,3

103:     test:
104:       suffix: tet
105:       requires: ctetgen

107:     test:
108:       suffix: tet_hash
109:       requires: ctetgen
110:       args: -dm_refine 1 -dm_plex_hash_location

112:     test:
113:       suffix: hex
114:       args: -dm_plex_simplex 0

116:     test:
117:       suffix: hex_hash
118:       args: -dm_plex_simplex 0 -dm_refine 1 -dm_plex_hash_location

120: TEST*/