1 static char help[] = "Tests for periodic mesh output\n\n"; 2 3 #include <petscdmplex.h> 4 5 typedef struct { 6 DM dm; 7 PetscInt dim; /* The topological mesh dimension */ 8 PetscBool cellSimplex; /* Use simplices or hexes */ 9 PetscInt faces[3]; /* Faces per direction */ 10 PetscBool isPeriodic; /* Flag for periodic mesh */ 11 DMBoundaryType periodicity[3]; /* Periodicity per direction */ 12 char filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */ 13 } AppCtx; 14 15 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 16 { 17 PetscInt n; 18 PetscErrorCode ierr; 19 20 PetscFunctionBegin; 21 options->dim = 2; 22 options->cellSimplex = PETSC_TRUE; 23 options->faces[0] = 1; 24 options->faces[1] = 1; 25 options->faces[2] = 1; 26 options->periodicity[0] = DM_BOUNDARY_NONE; 27 options->periodicity[1] = DM_BOUNDARY_NONE; 28 options->periodicity[2] = DM_BOUNDARY_NONE; 29 options->filename[0] = '\0'; 30 31 ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");CHKERRQ(ierr); 32 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex32.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr); 33 ierr = PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex32.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr); 34 n = 3; 35 ierr = PetscOptionsIntArray("-faces", "Faces per direction", "ex32.c", options->faces, &n, NULL);CHKERRQ(ierr); 36 n = 3; 37 ierr = PetscOptionsEnumArray("-periodicity", "Periodicity per direction", "ex32.c", DMBoundaryTypes, (PetscEnum *) options->periodicity, &n, &options->isPeriodic);CHKERRQ(ierr); 38 ierr = PetscOptionsString("-filename", "The mesh file", "ex32.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr); 39 ierr = PetscOptionsEnd(); 40 PetscFunctionReturn(0); 41 } 42 43 PetscErrorCode CheckMesh(DM dm, AppCtx *user) 44 { 45 PetscReal detJ, J[9], refVol = 1.0; 46 PetscReal vol; 47 PetscInt dim, depth, d, cStart, cEnd, c; 48 PetscErrorCode ierr; 49 50 PetscFunctionBegin; 51 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 52 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 53 for (d = 0; d < dim; ++d) { 54 refVol *= 2.0; 55 } 56 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 57 for (c = cStart; c < cEnd; ++c) { 58 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, NULL, J, NULL, &detJ);CHKERRQ(ierr); 59 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh cell %d is inverted, |J| = %g", c, detJ); 60 if (depth > 1) { 61 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr); 62 if (vol <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh cell %d is inverted, vol = %g", c, vol); 63 } 64 } 65 PetscFunctionReturn(0); 66 } 67 68 PetscErrorCode CompareCones(DM dm, DM idm) 69 { 70 PetscInt cStart, cEnd, c, vStart, vEnd, v; 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 75 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 76 for (c = cStart; c < cEnd; ++c) { 77 const PetscInt *cone; 78 PetscInt *points = NULL, numPoints, p, numVertices = 0, coneSize; 79 80 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 81 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 82 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 83 for (p = 0; p < numPoints*2; p += 2) { 84 const PetscInt point = points[p]; 85 if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point; 86 } 87 if (numVertices != coneSize) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In cell %d, cone size %d != %d vertices in closure", c, coneSize, numVertices); 88 for (v = 0; v < numVertices; ++v) { 89 if (cone[v] != points[v]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In cell %d, cone point %d is %d != %d vertex in closure", c, v, cone[v], points[v]); 90 } 91 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 92 } 93 PetscFunctionReturn(0); 94 } 95 96 PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 97 { 98 PetscInt dim = user->dim; 99 PetscBool cellSimplex = user->cellSimplex; 100 const char *filename = user->filename; 101 size_t len; 102 PetscMPIInt rank; 103 PetscErrorCode ierr; 104 105 PetscFunctionBegin; 106 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 107 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 108 if (len) { 109 ierr = DMPlexCreateFromFile(comm, filename, PETSC_FALSE, dm);CHKERRQ(ierr); 110 ierr = DMGetDimension(*dm, &dim);CHKERRQ(ierr); 111 } else { 112 PetscReal L[3] = {1.0, 1.0, 1.0}; 113 PetscReal maxCell[3]; 114 PetscInt d; 115 116 for (d = 0; d < dim; ++d) {maxCell[d] = (1.0/user->faces[d])*1.1;} 117 ierr = DMPlexCreateBoxMesh(comm, dim, cellSimplex, user->faces, NULL, NULL, user->periodicity, PETSC_TRUE, dm);CHKERRQ(ierr); 118 ierr = DMSetPeriodicity(*dm, user->isPeriodic, maxCell, L, user->periodicity);CHKERRQ(ierr); 119 } 120 { 121 DM pdm = NULL; 122 PetscPartitioner part; 123 124 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 125 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 126 /* Distribute mesh over processes */ 127 ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr); 128 if (pdm) { 129 ierr = DMDestroy(dm);CHKERRQ(ierr); 130 *dm = pdm; 131 } 132 } 133 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 134 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 135 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 136 user->dm = *dm; 137 PetscFunctionReturn(0); 138 } 139 140 int main(int argc, char **argv) 141 { 142 AppCtx user; /* user-defined work context */ 143 PetscErrorCode ierr; 144 145 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 146 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 147 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &user.dm);CHKERRQ(ierr); 148 ierr = CheckMesh(user.dm, &user);CHKERRQ(ierr); 149 ierr = DMDestroy(&user.dm);CHKERRQ(ierr); 150 ierr = PetscFinalize(); 151 return ierr; 152 } 153 154 /*TEST 155 156 test: 157 suffix: 0 158 args: -dim 2 -cell_simplex 0 -faces 3,1,0 -periodicity periodic,none,none -dm_view ::ascii_info_detail 159 test: 160 suffix: 1 161 nsize: 2 162 args: -dim 2 -cell_simplex 0 -faces 3,1,0 -periodicity periodic,none,none -petscpartitioner_type simple -dm_view ::ascii_info_detail 163 test: 164 suffix: 2 165 nsize: 2 166 args: -dim 2 -cell_simplex 0 -faces 6,2,0 -periodicity periodic,none,none -petscpartitioner_type simple -dm_view ::ascii_info_detail 167 test: 168 suffix: 3 169 nsize: 4 170 args: -dim 2 -cell_simplex 0 -faces 6,2,0 -periodicity periodic,none,none -petscpartitioner_type simple -dm_view ::ascii_info_detail 171 test: 172 suffix: 4 173 nsize: 2 174 args: -dim 2 -cell_simplex 0 -faces 3,1,0 -periodicity periodic,none,none -dm_plex_periodic_cut -petscpartitioner_type simple -dm_view ::ascii_info_detail 175 test: 176 suffix: 5 177 nsize: 2 178 args: -dim 2 -cell_simplex 0 -faces 6,2,0 -periodicity periodic,none,none -dm_plex_periodic_cut -petscpartitioner_type simple -dm_view ::ascii_info_detail 179 test: 180 suffix: 6 181 nsize: 4 182 args: -dim 2 -cell_simplex 0 -faces 6,2,0 -periodicity periodic,none,none -dm_plex_periodic_cut -petscpartitioner_type simple -dm_view ::ascii_info_detail 183 184 TEST*/ 185