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