1 static char help[] = "Tests interpolation and output of hybrid meshes\n\n"; 2 3 #include <petscdmplex.h> 4 5 typedef struct { 6 char filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */ 7 PetscBool interpolate; /* Interpolate the mesh */ 8 PetscInt meshNum; /* Which mesh we should construct */ 9 } AppCtx; 10 11 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 12 { 13 PetscErrorCode ierr; 14 15 PetscFunctionBeginUser; 16 options->filename[0] = '\0'; 17 options->interpolate = PETSC_FALSE; 18 options->meshNum = 0; 19 20 ierr = PetscOptionsBegin(comm, "", "Hybrid Output Test Options", "DMPLEX");CHKERRQ(ierr); 21 ierr = PetscOptionsString("-filename", "The mesh file", "ex8.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr); 22 ierr = PetscOptionsBool("-interpolate", "Interpolate the mesh", "ex8.c", options->interpolate, &options->interpolate, NULL);CHKERRQ(ierr); 23 ierr = PetscOptionsBoundedInt("-mesh_num", "The mesh we should construct", "ex8.c", options->meshNum, &options->meshNum, NULL,0);CHKERRQ(ierr); 24 ierr = PetscOptionsEnd(); 25 26 PetscFunctionReturn(0); 27 } 28 29 static PetscErrorCode CreateHybridMesh(MPI_Comm comm, PetscBool interpolate, DM *dm) 30 { 31 PetscInt dim; 32 PetscErrorCode ierr; 33 34 PetscFunctionBegin; 35 dim = 3; 36 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 37 ierr = PetscObjectSetName((PetscObject) *dm, "Simple Hybrid Mesh");CHKERRQ(ierr); 38 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 39 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 40 { 41 /* Simple mesh with 2 tets and 1 wedge */ 42 PetscInt numPoints[2] = {8, 3}; 43 PetscInt coneSize[11] = {4, 4, 6, 0, 0, 0, 0, 0, 0, 0, 0}; 44 PetscInt cones[14] = {4, 5, 6, 3, 7, 9, 8, 10, 4, 5, 6, 7, 8, 9}; 45 PetscInt coneOrientations[14] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 46 PetscScalar vertexCoords[48] = {-1.0, 1.0, 0.0, 47 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 0.0, 1.0, 1.0, 48 1.0, 0.0, 0.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 49 2.0, 1.0, 0.0}; 50 51 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 52 ierr = DMPlexSetHybridBounds(*dm, 2, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 53 if (interpolate) { 54 DM idm; 55 56 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 57 ierr = DMDestroy(dm);CHKERRQ(ierr); 58 *dm = idm; 59 } 60 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 61 } 62 PetscFunctionReturn(0); 63 } 64 65 /* TODO Need to extend interpolation to work when regular cells join to hybrid faces, rather than end faces */ 66 static PetscErrorCode CreateReverseHybridMesh(MPI_Comm comm, PetscBool interpolate, DM *dm) 67 { 68 PetscInt dim; 69 PetscErrorCode ierr; 70 71 PetscFunctionBegin; 72 dim = 3; 73 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 74 ierr = PetscObjectSetName((PetscObject) *dm, "Reverse Hybrid Mesh");CHKERRQ(ierr); 75 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 76 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 77 { 78 /* Simple mesh with 2 hexes and 3 wedges */ 79 PetscInt numPoints[2] = {16, 5}; 80 PetscInt coneSize[21] = {8, 8, 6, 6, 6, 81 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 82 PetscInt cones[34] = { 5, 6, 12, 11, 8, 14, 15, 9, 83 6, 7, 13, 12, 9, 15, 16, 10, 84 11, 17, 12, 14, 19, 15, 85 12, 18, 13, 15, 20, 16, 86 12, 17, 18, 15, 19, 20}; 87 PetscInt coneOrientations[34] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 88 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 89 PetscScalar vertexCoords[48] = {-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 1.0, 0.0, 90 -1.0, -1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 1.0, 1.0, 91 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 92 0.0, -1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 93 1.0, -1.0, 0.0, 1.0, 1.0, 0.0, 94 1.0, -1.0, 1.0, 1.0, 1.0, 1.0}; 95 96 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 97 ierr = DMPlexSetHybridBounds(*dm, 2, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 98 if (interpolate) { 99 DM idm; 100 101 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 102 ierr = DMDestroy(dm);CHKERRQ(ierr); 103 *dm = idm; 104 } 105 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 106 } 107 PetscFunctionReturn(0); 108 } 109 110 static PetscErrorCode OrderHybridMesh(DM *dm) 111 { 112 DM pdm; 113 IS perm; 114 PetscInt *ind; 115 PetscInt dim, pStart, pEnd, p, cStart, cEnd, c, Nhyb = 0, off[2]; 116 PetscErrorCode ierr; 117 118 PetscFunctionBegin; 119 ierr = DMGetDimension(*dm, &dim);CHKERRQ(ierr); 120 if (dim != 3) SETERRQ1(PetscObjectComm((PetscObject) *dm), PETSC_ERR_SUP, "No support for dimension %D", dim); 121 ierr = DMPlexGetChart(*dm, &pStart, &pEnd);CHKERRQ(ierr); 122 ierr = PetscMalloc1(pEnd-pStart, &ind);CHKERRQ(ierr); 123 for (p = 0; p < pEnd-pStart; ++p) ind[p] = p; 124 ierr = DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 125 for (c = cStart; c < cEnd; ++c) { 126 PetscInt coneSize; 127 128 ierr = DMPlexGetConeSize(*dm, c, &coneSize);CHKERRQ(ierr); 129 if (coneSize == 6) ++Nhyb; 130 } 131 off[0] = 0; 132 off[1] = cEnd - Nhyb; 133 for (c = cStart; c < cEnd; ++c) { 134 PetscInt coneSize; 135 136 ierr = DMPlexGetConeSize(*dm, c, &coneSize);CHKERRQ(ierr); 137 if (coneSize == 6) ind[c] = off[1]++; 138 else ind[c] = off[0]++; 139 } 140 if (off[0] != cEnd - Nhyb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of normal cells %D should be %D", off[0], cEnd - Nhyb); 141 if (off[1] != cEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of hybrid cells %D should be %D", off[1] - off[0], Nhyb); 142 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, ind, PETSC_OWN_POINTER, &perm);CHKERRQ(ierr); 143 ierr = DMPlexPermute(*dm, perm, &pdm);CHKERRQ(ierr); 144 ierr = DMPlexSetHybridBounds(pdm, cEnd-Nhyb, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 145 ierr = ISDestroy(&perm);CHKERRQ(ierr); 146 ierr = DMDestroy(dm);CHKERRQ(ierr); 147 *dm = pdm; 148 PetscFunctionReturn(0); 149 } 150 151 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 152 { 153 const char *filename = user->filename; 154 PetscBool interpolate = user->interpolate; 155 PetscInt meshNum = user->meshNum; 156 size_t len; 157 PetscErrorCode ierr; 158 159 PetscFunctionBegin; 160 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 161 if (len) { 162 ierr = DMPlexCreateFromFile(comm, filename, PETSC_FALSE, dm);CHKERRQ(ierr); 163 ierr = OrderHybridMesh(dm);CHKERRQ(ierr); 164 if (interpolate) { 165 DM idm; 166 167 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 168 ierr = DMDestroy(dm);CHKERRQ(ierr); 169 *dm = idm; 170 } 171 ierr = PetscObjectSetName((PetscObject) *dm, "Input Mesh");CHKERRQ(ierr); 172 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 173 } else { 174 switch (meshNum) { 175 case 0: 176 ierr = CreateHybridMesh(comm, interpolate, dm);CHKERRQ(ierr);break; 177 case 1: 178 ierr = CreateReverseHybridMesh(comm, interpolate, dm);CHKERRQ(ierr);break; 179 default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Unknown mesh number %D", user->meshNum); 180 } 181 } 182 PetscFunctionReturn(0); 183 } 184 185 int main(int argc, char **argv) 186 { 187 DM dm; 188 AppCtx user; 189 PetscErrorCode ierr; 190 191 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 192 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 193 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 194 ierr = DMDestroy(&dm);CHKERRQ(ierr); 195 ierr = PetscFinalize(); 196 return ierr; 197 } 198 199 /*TEST 200 201 test: 202 suffix: 0 203 args: -interpolate -dm_view ascii::ascii_info_detail 204 205 test: 206 suffix: 1 207 args: -mesh_num 1 -interpolate -dm_view ascii::ascii_info_detail 208 209 TEST*/ 210