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