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