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