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