Actual source code: ex20.c

      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";

      6: #include <petsc.h>
      7: #include <petscdm.h>
      8: #include <petscdmda.h>
      9: #include <petscdmplex.h>
     10: #include <petscdmswarm.h>

     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];

     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));

     26:   PetscCall(DMDASetElementType(celldm, DMDA_ELEMENT_Q1));
     27:   PetscCall(DMSetFromOptions(celldm));
     28:   PetscCall(DMSetUp(celldm));

     30:   PetscCall(DMDASetUniformCoordinates(celldm, 0.0, 2.0, 0.0, 1.0, 0.0, 1.5));

     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));

     38:   /* Configure swarm to be of type PIC */
     39:   PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
     40:   PetscCall(DMSwarmSetCellDM(swarm, celldm));

     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));

     47:   /* Set initial local sizes of the DMSwarm with a buffer length of zero */
     48:   PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));

     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));

     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));

     66:   {
     67:     PetscInt    npoints, *list;
     68:     PetscMPIInt rank;

     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: }

     83: PetscErrorCode pic_insert_DMPLEX_with_cell_list(PetscInt dim)
     84: {
     85:   DM          celldm = NULL, swarm, distributedMesh = NULL;
     86:   const char *fieldnames[] = {"viscosity"};

     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;

     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];

    105:     PetscCall(PetscMalloc1(n_tricells * 3, &tricells));
    106:     PetscCall(PetscMalloc1(nx[0] * nx[1] * 2, &trivert));

    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:     }

    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;

    126:         idx  = (ii) + (jj)*nx[0];
    127:         idx0 = idx;
    128:         idx1 = idx0 + 1;
    129:         idx2 = idx1 + nx[0];
    130:         idx3 = idx0 + nx[0];

    132:         tricells[3 * cnt + 0] = idx0;
    133:         tricells[3 * cnt + 1] = idx1;
    134:         tricells[3 * cnt + 2] = idx2;
    135:         cnt++;

    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");

    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;

    163:     PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &section));
    164:     PetscCall(DMSetLocalSection(celldm, section));
    165:     PetscCall(PetscSectionDestroy(&section));
    166:   }
    167:   PetscCall(DMSetUp(celldm));
    168:   {
    169:     PetscViewer viewer;

    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:   }

    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));

    185:   PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
    186:   PetscCall(DMSwarmSetCellDM(swarm, celldm));

    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));

    193:   /* Set initial local sizes of the DMSwarm with a buffer length of zero */
    194:   PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));

    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: }

    206: PetscErrorCode pic_insert_DMPLEX(PetscBool is_simplex, PetscInt dim)
    207: {
    208:   DM          celldm, swarm, distributedMesh = NULL;
    209:   const char *fieldnames[] = {"viscosity", "DMSwarm_rank"};

    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, 0, PETSC_TRUE, &celldm));
    216:   }

    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;

    232:     PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &section));
    233:     PetscCall(DMSetLocalSection(celldm, section));
    234:     PetscCall(PetscSectionDestroy(&section));
    235:   }
    236:   PetscCall(DMSetUp(celldm));
    237:   {
    238:     PetscViewer viewer;

    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:   }

    248:   PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm));
    249:   PetscCall(PetscObjectSetName((PetscObject)swarm, "Swarm"));
    250:   PetscCall(DMSetType(swarm, DMSWARM));
    251:   PetscCall(DMSetDimension(swarm, dim));

    253:   PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
    254:   PetscCall(DMSwarmSetCellDM(swarm, celldm));

    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));

    261:   /* Set initial local sizes of the DMSwarm with a buffer length of zero */
    262:   PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));

    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: }

    274: int main(int argc, char **args)
    275: {
    276:   PetscInt mode = 0;
    277:   PetscInt dim  = 2;

    279:   PetscFunctionBeginUser;
    280:   PetscCall(PetscInitialize(&argc, &args, NULL, 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: }

    303: /*TEST

    305:    test:
    306:       args:
    307:       requires: !complex double
    308:       filter: grep -v atomic
    309:       filter_output: grep -v atomic

    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

    318: TEST*/