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