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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMDASetElementType(celldm,DMDA_ELEMENT_Q1)); 32 CHKERRQ(DMSetFromOptions(celldm)); 33 CHKERRQ(DMSetUp(celldm)); 34 35 CHKERRQ(DMDASetUniformCoordinates(celldm,0.0,2.0,0.0,1.0,0.0,1.5)); 36 37 /* Create the DMSwarm */ 38 CHKERRQ(DMCreate(PETSC_COMM_WORLD,&swarm)); 39 CHKERRQ(PetscObjectSetName((PetscObject)swarm,"Swarm")); 40 CHKERRQ(DMSetType(swarm,DMSWARM)); 41 CHKERRQ(DMSetDimension(swarm,dim)); 42 43 /* Configure swarm to be of type PIC */ 44 CHKERRQ(DMSwarmSetType(swarm,DMSWARM_PIC)); 45 CHKERRQ(DMSwarmSetCellDM(swarm,celldm)); 46 47 /* Register two scalar fields within the DMSwarm */ 48 CHKERRQ(DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE)); 49 CHKERRQ(DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE)); 50 CHKERRQ(DMSwarmFinalizeFieldRegister(swarm)); 51 52 /* Set initial local sizes of the DMSwarm with a buffer length of zero */ 53 CHKERRQ(DMSwarmSetLocalSizes(swarm,4,0)); 54 55 /* Insert swarm coordinates cell-wise */ 56 CHKERRQ(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 CHKERRQ(DMSwarmSetPointsUniformCoordinates(swarm,min,max,ndir,ADD_VALUES)); 62 63 /* This should be dispatched from a regular DMView() */ 64 CHKERRQ(DMSwarmViewXDMF(swarm,"ex20.xmf")); 65 CHKERRQ(DMView(celldm,PETSC_VIEWER_STDOUT_WORLD)); 66 CHKERRQ(DMView(swarm,PETSC_VIEWER_STDOUT_WORLD)); 67 68 { 69 PetscInt npoints,*list; 70 PetscMPIInt rank; 71 72 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 73 CHKERRQ(DMSwarmSortGetAccess(swarm)); 74 CHKERRQ(DMSwarmSortGetNumberOfPointsPerCell(swarm,0,&npoints)); 75 CHKERRQ(DMSwarmSortGetPointsPerCell(swarm,rank,&npoints,&list)); 76 CHKERRQ(PetscFree(list)); 77 CHKERRQ(DMSwarmSortRestoreAccess(swarm)); 78 } 79 CHKERRQ(DMSwarmMigrate(swarm,PETSC_FALSE)); 80 CHKERRQ(DMDestroy(&celldm)); 81 CHKERRQ(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 CHKERRQ(PetscMalloc1(n_tricells*3,&tricells)); 108 CHKERRQ(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 CHKERRQ(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD,dim,n_tricells,n_trivert,3,PETSC_TRUE,tricells,dim,trivert,&celldm)); 146 CHKERRQ(PetscFree(trivert)); 147 CHKERRQ(PetscFree(tricells)); 148 } 149 PetscCheckFalse(dim == 3,PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only 2D PLEX example supported"); 150 151 /* Distribute mesh over processes */ 152 CHKERRQ(DMPlexDistribute(celldm,0,NULL,&distributedMesh)); 153 if (distributedMesh) { 154 CHKERRQ(DMDestroy(&celldm)); 155 celldm = distributedMesh; 156 } 157 CHKERRQ(PetscObjectSetName((PetscObject)celldm,"Cells")); 158 CHKERRQ(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 CHKERRQ(DMPlexCreateSection(celldm,NULL,numComp,numDof,numBC,NULL,NULL,NULL,NULL,§ion)); 166 CHKERRQ(DMSetLocalSection(celldm,section)); 167 CHKERRQ(PetscSectionDestroy(§ion)); 168 } 169 CHKERRQ(DMSetUp(celldm)); 170 { 171 PetscViewer viewer; 172 173 CHKERRQ(PetscViewerCreate(PETSC_COMM_WORLD,&viewer)); 174 CHKERRQ(PetscViewerSetType(viewer,PETSCVIEWERVTK)); 175 CHKERRQ(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE)); 176 CHKERRQ(PetscViewerFileSetName(viewer,"ex20plex.vtk")); 177 CHKERRQ(DMView(celldm,viewer)); 178 CHKERRQ(PetscViewerDestroy(&viewer)); 179 } 180 181 /* Create the DMSwarm */ 182 CHKERRQ(DMCreate(PETSC_COMM_WORLD,&swarm)); 183 CHKERRQ(PetscObjectSetName((PetscObject)swarm,"Swarm")); 184 CHKERRQ(DMSetType(swarm,DMSWARM)); 185 CHKERRQ(DMSetDimension(swarm,dim)); 186 187 CHKERRQ(DMSwarmSetType(swarm,DMSWARM_PIC)); 188 CHKERRQ(DMSwarmSetCellDM(swarm,celldm)); 189 190 /* Register two scalar fields within the DMSwarm */ 191 CHKERRQ(DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE)); 192 CHKERRQ(DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE)); 193 CHKERRQ(DMSwarmFinalizeFieldRegister(swarm)); 194 195 /* Set initial local sizes of the DMSwarm with a buffer length of zero */ 196 CHKERRQ(DMSwarmSetLocalSizes(swarm,4,0)); 197 198 /* Insert swarm coordinates cell-wise */ 199 CHKERRQ(DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_SUBDIVISION,2)); 200 CHKERRQ(DMSwarmViewFieldsXDMF(swarm,"ex20.xmf",1,fieldnames)); 201 CHKERRQ(DMView(celldm,PETSC_VIEWER_STDOUT_WORLD)); 202 CHKERRQ(DMView(swarm,PETSC_VIEWER_STDOUT_WORLD)); 203 CHKERRQ(DMDestroy(&celldm)); 204 CHKERRQ(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 CHKERRQ(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, is_simplex, faces, NULL, NULL, NULL, PETSC_TRUE, &celldm)); 219 } 220 221 /* Distribute mesh over processes */ 222 CHKERRQ(DMPlexDistribute(celldm,0,NULL,&distributedMesh)); 223 if (distributedMesh) { 224 CHKERRQ(DMDestroy(&celldm)); 225 celldm = distributedMesh; 226 } 227 CHKERRQ(PetscObjectSetName((PetscObject)celldm,"Cells")); 228 CHKERRQ(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 CHKERRQ(DMPlexCreateSection(celldm,NULL,numComp,numDof,numBC,NULL,NULL,NULL,NULL,§ion)); 236 CHKERRQ(DMSetLocalSection(celldm,section)); 237 CHKERRQ(PetscSectionDestroy(§ion)); 238 } 239 CHKERRQ(DMSetUp(celldm)); 240 { 241 PetscViewer viewer; 242 243 CHKERRQ(PetscViewerCreate(PETSC_COMM_WORLD,&viewer)); 244 CHKERRQ(PetscViewerSetType(viewer,PETSCVIEWERVTK)); 245 CHKERRQ(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE)); 246 CHKERRQ(PetscViewerFileSetName(viewer,"ex20plex.vtk")); 247 CHKERRQ(DMView(celldm,viewer)); 248 CHKERRQ(PetscViewerDestroy(&viewer)); 249 } 250 251 CHKERRQ(DMCreate(PETSC_COMM_WORLD,&swarm)); 252 CHKERRQ(PetscObjectSetName((PetscObject)swarm,"Swarm")); 253 CHKERRQ(DMSetType(swarm,DMSWARM)); 254 CHKERRQ(DMSetDimension(swarm,dim)); 255 256 CHKERRQ(DMSwarmSetType(swarm,DMSWARM_PIC)); 257 CHKERRQ(DMSwarmSetCellDM(swarm,celldm)); 258 259 /* Register two scalar fields within the DMSwarm */ 260 CHKERRQ(DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE)); 261 CHKERRQ(DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE)); 262 CHKERRQ(DMSwarmFinalizeFieldRegister(swarm)); 263 264 /* Set initial local sizes of the DMSwarm with a buffer length of zero */ 265 CHKERRQ(DMSwarmSetLocalSizes(swarm,4,0)); 266 267 /* Insert swarm coordinates cell-wise */ 268 CHKERRQ(DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_GAUSS,3)); 269 CHKERRQ(DMSwarmViewFieldsXDMF(swarm,"ex20.xmf",2,fieldnames)); 270 CHKERRQ(DMView(celldm,PETSC_VIEWER_STDOUT_WORLD)); 271 CHKERRQ(DMView(swarm,PETSC_VIEWER_STDOUT_WORLD)); 272 CHKERRQ(DMDestroy(&celldm)); 273 CHKERRQ(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 CHKERRQ(PetscInitialize(&argc,&args,(char*)0,help)); 283 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mode",&mode,NULL)); 284 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL)); 285 switch (mode) { 286 case 0: 287 CHKERRQ(pic_insert_DMDA(dim)); 288 break; 289 case 1: 290 /* tri / tet */ 291 CHKERRQ(pic_insert_DMPLEX(PETSC_TRUE,dim)); 292 break; 293 case 2: 294 /* quad / hex */ 295 CHKERRQ(pic_insert_DMPLEX(PETSC_FALSE,dim)); 296 break; 297 default: 298 CHKERRQ(pic_insert_DMDA(dim)); 299 break; 300 } 301 CHKERRQ(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