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