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