1 2 static char help[] = "DMSwarm-PIC demonstator of inserting points into a cell DM \n\ 3 Options: \n\ 4 -mode {0,1} : 0 ==> DMDA, 1 ==> DMPLEX cell DM \n\ 5 -dim {2,3} : spatial dimension\n"; 6 7 #include <petsc.h> 8 #include <petscdm.h> 9 #include <petscdmda.h> 10 #include <petscdmplex.h> 11 #include <petscdmswarm.h> 12 13 PetscErrorCode pic_insert_DMDA(PetscInt dim) 14 { 15 DM celldm = NULL, swarm; 16 PetscInt dof, stencil_width; 17 PetscReal min[3], max[3]; 18 PetscInt ndir[3]; 19 20 PetscFunctionBegin; 21 /* Create the background cell DM */ 22 dof = 1; 23 stencil_width = 1; 24 if (dim == 2) PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 25, 13, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, &celldm)); 25 if (dim == 3) PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 25, 13, 19, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, NULL, &celldm)); 26 27 PetscCall(DMDASetElementType(celldm, DMDA_ELEMENT_Q1)); 28 PetscCall(DMSetFromOptions(celldm)); 29 PetscCall(DMSetUp(celldm)); 30 31 PetscCall(DMDASetUniformCoordinates(celldm, 0.0, 2.0, 0.0, 1.0, 0.0, 1.5)); 32 33 /* Create the DMSwarm */ 34 PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm)); 35 PetscCall(PetscObjectSetName((PetscObject)swarm, "Swarm")); 36 PetscCall(DMSetType(swarm, DMSWARM)); 37 PetscCall(DMSetDimension(swarm, dim)); 38 39 /* Configure swarm to be of type PIC */ 40 PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC)); 41 PetscCall(DMSwarmSetCellDM(swarm, celldm)); 42 43 /* Register two scalar fields within the DMSwarm */ 44 PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE)); 45 PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE)); 46 PetscCall(DMSwarmFinalizeFieldRegister(swarm)); 47 48 /* Set initial local sizes of the DMSwarm with a buffer length of zero */ 49 PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0)); 50 51 /* Insert swarm coordinates cell-wise */ 52 PetscCall(DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_REGULAR, 3)); 53 min[0] = 0.5; 54 max[0] = 0.7; 55 min[1] = 0.5; 56 max[1] = 0.8; 57 min[2] = 0.5; 58 max[2] = 0.9; 59 ndir[0] = ndir[1] = ndir[2] = 30; 60 PetscCall(DMSwarmSetPointsUniformCoordinates(swarm, min, max, ndir, ADD_VALUES)); 61 62 /* This should be dispatched from a regular DMView() */ 63 PetscCall(DMSwarmViewXDMF(swarm, "ex20.xmf")); 64 PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD)); 65 PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD)); 66 67 { 68 PetscInt npoints, *list; 69 PetscMPIInt rank; 70 71 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 72 PetscCall(DMSwarmSortGetAccess(swarm)); 73 PetscCall(DMSwarmSortGetNumberOfPointsPerCell(swarm, 0, &npoints)); 74 PetscCall(DMSwarmSortGetPointsPerCell(swarm, rank, &npoints, &list)); 75 PetscCall(PetscFree(list)); 76 PetscCall(DMSwarmSortRestoreAccess(swarm)); 77 } 78 PetscCall(DMSwarmMigrate(swarm, PETSC_FALSE)); 79 PetscCall(DMDestroy(&celldm)); 80 PetscCall(DMDestroy(&swarm)); 81 PetscFunctionReturn(0); 82 } 83 84 PetscErrorCode pic_insert_DMPLEX_with_cell_list(PetscInt dim) 85 { 86 DM celldm = NULL, swarm, distributedMesh = NULL; 87 const char *fieldnames[] = {"viscosity"}; 88 89 PetscFunctionBegin; 90 /* Create the background cell DM */ 91 if (dim == 2) { 92 PetscInt cells_per_dim[2], nx[2]; 93 PetscInt n_tricells; 94 PetscInt n_trivert; 95 PetscInt *tricells; 96 PetscReal *trivert, dx, dy; 97 PetscInt ii, jj, cnt; 98 99 cells_per_dim[0] = 4; 100 cells_per_dim[1] = 4; 101 n_tricells = cells_per_dim[0] * cells_per_dim[1] * 2; 102 nx[0] = cells_per_dim[0] + 1; 103 nx[1] = cells_per_dim[1] + 1; 104 n_trivert = nx[0] * nx[1]; 105 106 PetscCall(PetscMalloc1(n_tricells * 3, &tricells)); 107 PetscCall(PetscMalloc1(nx[0] * nx[1] * 2, &trivert)); 108 109 /* verts */ 110 cnt = 0; 111 dx = 2.0 / ((PetscReal)cells_per_dim[0]); 112 dy = 1.0 / ((PetscReal)cells_per_dim[1]); 113 for (jj = 0; jj < nx[1]; jj++) { 114 for (ii = 0; ii < nx[0]; ii++) { 115 trivert[2 * cnt + 0] = 0.0 + ii * dx; 116 trivert[2 * cnt + 1] = 0.0 + jj * dy; 117 cnt++; 118 } 119 } 120 121 /* connectivity */ 122 cnt = 0; 123 for (jj = 0; jj < cells_per_dim[1]; jj++) { 124 for (ii = 0; ii < cells_per_dim[0]; ii++) { 125 PetscInt idx, idx0, idx1, idx2, idx3; 126 127 idx = (ii) + (jj)*nx[0]; 128 idx0 = idx; 129 idx1 = idx0 + 1; 130 idx2 = idx1 + nx[0]; 131 idx3 = idx0 + nx[0]; 132 133 tricells[3 * cnt + 0] = idx0; 134 tricells[3 * cnt + 1] = idx1; 135 tricells[3 * cnt + 2] = idx2; 136 cnt++; 137 138 tricells[3 * cnt + 0] = idx0; 139 tricells[3 * cnt + 1] = idx2; 140 tricells[3 * cnt + 2] = idx3; 141 cnt++; 142 } 143 } 144 PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, n_tricells, n_trivert, 3, PETSC_TRUE, tricells, dim, trivert, &celldm)); 145 PetscCall(PetscFree(trivert)); 146 PetscCall(PetscFree(tricells)); 147 } 148 PetscCheck(dim != 3, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Only 2D PLEX example supported"); 149 150 /* Distribute mesh over processes */ 151 PetscCall(DMPlexDistribute(celldm, 0, NULL, &distributedMesh)); 152 if (distributedMesh) { 153 PetscCall(DMDestroy(&celldm)); 154 celldm = distributedMesh; 155 } 156 PetscCall(PetscObjectSetName((PetscObject)celldm, "Cells")); 157 PetscCall(DMSetFromOptions(celldm)); 158 { 159 PetscInt numComp[] = {1}; 160 PetscInt numDof[] = {1, 0, 0}; /* vert, edge, cell */ 161 PetscInt numBC = 0; 162 PetscSection section; 163 164 PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, §ion)); 165 PetscCall(DMSetLocalSection(celldm, section)); 166 PetscCall(PetscSectionDestroy(§ion)); 167 } 168 PetscCall(DMSetUp(celldm)); 169 { 170 PetscViewer viewer; 171 172 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer)); 173 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK)); 174 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE)); 175 PetscCall(PetscViewerFileSetName(viewer, "ex20plex.vtk")); 176 PetscCall(DMView(celldm, viewer)); 177 PetscCall(PetscViewerDestroy(&viewer)); 178 } 179 180 /* Create the DMSwarm */ 181 PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm)); 182 PetscCall(PetscObjectSetName((PetscObject)swarm, "Swarm")); 183 PetscCall(DMSetType(swarm, DMSWARM)); 184 PetscCall(DMSetDimension(swarm, dim)); 185 186 PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC)); 187 PetscCall(DMSwarmSetCellDM(swarm, celldm)); 188 189 /* Register two scalar fields within the DMSwarm */ 190 PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE)); 191 PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE)); 192 PetscCall(DMSwarmFinalizeFieldRegister(swarm)); 193 194 /* Set initial local sizes of the DMSwarm with a buffer length of zero */ 195 PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0)); 196 197 /* Insert swarm coordinates cell-wise */ 198 PetscCall(DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_SUBDIVISION, 2)); 199 PetscCall(DMSwarmViewFieldsXDMF(swarm, "ex20.xmf", 1, fieldnames)); 200 PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD)); 201 PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD)); 202 PetscCall(DMDestroy(&celldm)); 203 PetscCall(DMDestroy(&swarm)); 204 PetscFunctionReturn(0); 205 } 206 207 PetscErrorCode pic_insert_DMPLEX(PetscBool is_simplex, PetscInt dim) 208 { 209 DM celldm, swarm, distributedMesh = NULL; 210 const char *fieldnames[] = {"viscosity", "DMSwarm_rank"}; 211 212 PetscFunctionBegin; 213 214 /* Create the background cell DM */ 215 { 216 PetscInt faces[3] = {4, 2, 4}; 217 PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, is_simplex, faces, NULL, NULL, NULL, PETSC_TRUE, &celldm)); 218 } 219 220 /* Distribute mesh over processes */ 221 PetscCall(DMPlexDistribute(celldm, 0, NULL, &distributedMesh)); 222 if (distributedMesh) { 223 PetscCall(DMDestroy(&celldm)); 224 celldm = distributedMesh; 225 } 226 PetscCall(PetscObjectSetName((PetscObject)celldm, "Cells")); 227 PetscCall(DMSetFromOptions(celldm)); 228 { 229 PetscInt numComp[] = {1}; 230 PetscInt numDof[] = {1, 0, 0}; /* vert, edge, cell */ 231 PetscInt numBC = 0; 232 PetscSection section; 233 234 PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, §ion)); 235 PetscCall(DMSetLocalSection(celldm, section)); 236 PetscCall(PetscSectionDestroy(§ion)); 237 } 238 PetscCall(DMSetUp(celldm)); 239 { 240 PetscViewer viewer; 241 242 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer)); 243 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK)); 244 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE)); 245 PetscCall(PetscViewerFileSetName(viewer, "ex20plex.vtk")); 246 PetscCall(DMView(celldm, viewer)); 247 PetscCall(PetscViewerDestroy(&viewer)); 248 } 249 250 PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm)); 251 PetscCall(PetscObjectSetName((PetscObject)swarm, "Swarm")); 252 PetscCall(DMSetType(swarm, DMSWARM)); 253 PetscCall(DMSetDimension(swarm, dim)); 254 255 PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC)); 256 PetscCall(DMSwarmSetCellDM(swarm, celldm)); 257 258 /* Register two scalar fields within the DMSwarm */ 259 PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE)); 260 PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE)); 261 PetscCall(DMSwarmFinalizeFieldRegister(swarm)); 262 263 /* Set initial local sizes of the DMSwarm with a buffer length of zero */ 264 PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0)); 265 266 /* Insert swarm coordinates cell-wise */ 267 PetscCall(DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_GAUSS, 3)); 268 PetscCall(DMSwarmViewFieldsXDMF(swarm, "ex20.xmf", 2, fieldnames)); 269 PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD)); 270 PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD)); 271 PetscCall(DMDestroy(&celldm)); 272 PetscCall(DMDestroy(&swarm)); 273 PetscFunctionReturn(0); 274 } 275 276 int main(int argc, char **args) 277 { 278 PetscInt mode = 0; 279 PetscInt dim = 2; 280 281 PetscFunctionBeginUser; 282 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 283 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mode", &mode, NULL)); 284 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL)); 285 switch (mode) { 286 case 0: 287 PetscCall(pic_insert_DMDA(dim)); 288 break; 289 case 1: 290 /* tri / tet */ 291 PetscCall(pic_insert_DMPLEX(PETSC_TRUE, dim)); 292 break; 293 case 2: 294 /* quad / hex */ 295 PetscCall(pic_insert_DMPLEX(PETSC_FALSE, dim)); 296 break; 297 default: 298 PetscCall(pic_insert_DMDA(dim)); 299 break; 300 } 301 PetscCall(PetscFinalize()); 302 return 0; 303 } 304 305 /*TEST 306 307 test: 308 args: 309 requires: !complex double 310 filter: grep -v atomic 311 filter_output: grep -v atomic 312 313 test: 314 suffix: 2 315 requires: triangle double !complex 316 args: -mode 1 317 filter: grep -v atomic 318 filter_output: grep -v atomic 319 320 TEST*/ 321