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, sizeof(options->filename), 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();CHKERRQ(ierr); 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 if (interpolate) { 53 DM idm; 54 55 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 56 ierr = DMDestroy(dm);CHKERRQ(ierr); 57 *dm = idm; 58 } 59 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 60 } 61 PetscFunctionReturn(0); 62 } 63 64 /* 65 This is not a valid mesh. We need to either change to tensor quad prisms or regular triangular prisms. 66 67 10-------16--------20 68 /| | 69 / | | 70 / | | 71 9---|---15 | 72 /| 7 | 13--------18 73 / | / | / ____/ 74 / | / | /____/ 75 8 |/ 14---|//---19 76 | 6 | 12 77 | / | / \ 78 | / | / \__ 79 |/ |/ \ 80 5--------11--------17 81 */ 82 static PetscErrorCode CreateReverseHybridMesh(MPI_Comm comm, PetscBool interpolate, DM *dm) 83 { 84 PetscInt dim; 85 PetscErrorCode ierr; 86 87 PetscFunctionBegin; 88 dim = 3; 89 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 90 ierr = PetscObjectSetName((PetscObject) *dm, "Reverse Hybrid Mesh");CHKERRQ(ierr); 91 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 92 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 93 { 94 /* Simple mesh with 2 hexes and 3 wedges */ 95 PetscInt numPoints[2] = {16, 5}; 96 PetscInt coneSize[21] = {8, 8, 6, 6, 6, 97 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 98 PetscInt cones[34] = { 5, 6, 12, 11, 8, 14, 15, 9, 99 6, 7, 13, 12, 9, 15, 16, 10, 100 11, 17, 12, 14, 19, 15, 101 12, 18, 13, 15, 20, 16, 102 12, 17, 18, 15, 19, 20}; 103 PetscInt coneOrientations[34] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 104 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 105 PetscScalar vertexCoords[48] = {-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 1.0, 0.0, 106 -1.0, -1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 1.0, 1.0, 107 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 108 0.0, -1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 109 1.0, -1.0, 0.0, 1.0, 1.0, 0.0, 110 1.0, -1.0, 1.0, 1.0, 1.0, 1.0}; 111 112 ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 113 if (interpolate) { 114 DM idm; 115 116 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 117 ierr = DMDestroy(dm);CHKERRQ(ierr); 118 *dm = idm; 119 } 120 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 121 } 122 PetscFunctionReturn(0); 123 } 124 125 static PetscErrorCode OrderHybridMesh(DM *dm) 126 { 127 DM pdm; 128 IS perm; 129 PetscInt *ind; 130 PetscInt dim, pStart, pEnd, p, cStart, cEnd, c, Nhyb = 0, off[2]; 131 PetscErrorCode ierr; 132 133 PetscFunctionBegin; 134 ierr = DMGetDimension(*dm, &dim);CHKERRQ(ierr); 135 if (dim != 3) SETERRQ1(PetscObjectComm((PetscObject) *dm), PETSC_ERR_SUP, "No support for dimension %D", dim); 136 ierr = DMPlexGetChart(*dm, &pStart, &pEnd);CHKERRQ(ierr); 137 ierr = PetscMalloc1(pEnd-pStart, &ind);CHKERRQ(ierr); 138 for (p = 0; p < pEnd-pStart; ++p) ind[p] = p; 139 ierr = DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 140 for (c = cStart; c < cEnd; ++c) { 141 PetscInt coneSize; 142 143 ierr = DMPlexGetConeSize(*dm, c, &coneSize);CHKERRQ(ierr); 144 if (coneSize == 6) ++Nhyb; 145 } 146 off[0] = 0; 147 off[1] = cEnd - Nhyb; 148 for (c = cStart; c < cEnd; ++c) { 149 PetscInt coneSize; 150 151 ierr = DMPlexGetConeSize(*dm, c, &coneSize);CHKERRQ(ierr); 152 if (coneSize == 6) ind[c] = off[1]++; 153 else ind[c] = off[0]++; 154 } 155 if (off[0] != cEnd - Nhyb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of normal cells %D should be %D", off[0], cEnd - Nhyb); 156 if (off[1] != cEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of hybrid cells %D should be %D", off[1] - off[0], Nhyb); 157 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, ind, PETSC_OWN_POINTER, &perm);CHKERRQ(ierr); 158 ierr = DMPlexPermute(*dm, perm, &pdm);CHKERRQ(ierr); 159 ierr = ISDestroy(&perm);CHKERRQ(ierr); 160 ierr = DMDestroy(dm);CHKERRQ(ierr); 161 *dm = pdm; 162 PetscFunctionReturn(0); 163 } 164 165 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 166 { 167 const char *filename = user->filename; 168 PetscBool interpolate = user->interpolate; 169 PetscInt meshNum = user->meshNum; 170 size_t len; 171 PetscErrorCode ierr; 172 173 PetscFunctionBegin; 174 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 175 if (len) { 176 ierr = DMPlexCreateFromFile(comm, filename, PETSC_FALSE, dm);CHKERRQ(ierr); 177 ierr = OrderHybridMesh(dm);CHKERRQ(ierr); 178 if (interpolate) { 179 DM idm; 180 181 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 182 ierr = DMDestroy(dm);CHKERRQ(ierr); 183 *dm = idm; 184 } 185 ierr = PetscObjectSetName((PetscObject) *dm, "Input Mesh");CHKERRQ(ierr); 186 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 187 } else { 188 switch (meshNum) { 189 case 0: 190 ierr = CreateHybridMesh(comm, interpolate, dm);CHKERRQ(ierr);break; 191 case 1: 192 ierr = CreateReverseHybridMesh(comm, interpolate, dm);CHKERRQ(ierr);break; 193 default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Unknown mesh number %D", user->meshNum); 194 } 195 } 196 PetscFunctionReturn(0); 197 } 198 199 int main(int argc, char **argv) 200 { 201 DM dm; 202 AppCtx user; 203 PetscErrorCode ierr; 204 205 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 206 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 207 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 208 ierr = DMDestroy(&dm);CHKERRQ(ierr); 209 ierr = PetscFinalize(); 210 return ierr; 211 } 212 213 /*TEST 214 215 test: 216 suffix: 0 217 args: -interpolate -dm_view ascii::ascii_info_detail 218 219 # Test needs to be reworked 220 test: 221 requires: BROKEN 222 suffix: 1 223 args: -mesh_num 1 -interpolate -dm_view ascii::ascii_info_detail 224 225 TEST*/ 226