Actual source code: ex8.c

      1: static char help[] = "Element closure restrictions in tensor/lexicographic/spectral-element ordering using DMPlex\n\n";

      3: #include <petscdmplex.h>
      4: #include <petscds.h>

      6: static PetscErrorCode ViewOffsets(DM dm, Vec X)
      7: {
      8:   PetscInt           num_elem, elem_size, num_comp, num_dof;
      9:   PetscInt          *elem_restr_offsets;
     10:   const PetscScalar *x = NULL;
     11:   const char        *name;

     13:   PetscFunctionBegin;
     14:   PetscCall(PetscObjectGetName((PetscObject)dm, &name));
     15:   PetscCall(DMPlexGetLocalOffsets(dm, NULL, 0, 0, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets));
     16:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "DM %s offsets: num_elem %" PetscInt_FMT ", size %" PetscInt_FMT ", comp %" PetscInt_FMT ", dof %" PetscInt_FMT "\n", name, num_elem, elem_size, num_comp, num_dof));
     17:   if (X) PetscCall(VecGetArrayRead(X, &x));
     18:   for (PetscInt c = 0; c < num_elem; c++) {
     19:     PetscCall(PetscIntView(elem_size, &elem_restr_offsets[c * elem_size], PETSC_VIEWER_STDOUT_SELF));
     20:     if (x) {
     21:       for (PetscInt i = 0; i < elem_size; i++) PetscCall(PetscScalarView(num_comp, &x[elem_restr_offsets[c * elem_size + i]], PETSC_VIEWER_STDERR_SELF));
     22:     }
     23:   }
     24:   if (X) PetscCall(VecRestoreArrayRead(X, &x));
     25:   PetscCall(PetscFree(elem_restr_offsets));
     26:   PetscFunctionReturn(PETSC_SUCCESS);
     27: }

     29: int main(int argc, char **argv)
     30: {
     31:   DM           dm;
     32:   PetscSection section;
     33:   PetscFE      fe;
     34:   PetscInt     dim, Nf = 1, c, cStart, cEnd;
     35:   PetscBool    view_coord = PETSC_FALSE, tensor = PETSC_TRUE, project = PETSC_FALSE;

     37:   PetscFunctionBeginUser;
     38:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     39:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Tensor closure restrictions", "DMPLEX");
     40:   PetscCall(PetscOptionsBool("-closure_tensor", "Apply DMPlexSetClosurePermutationTensor", "ex8.c", tensor, &tensor, NULL));
     41:   PetscCall(PetscOptionsBool("-project_coordinates", "Call DMSetCoordinateDisc() explicitly", "ex8.c", project, &project, NULL));
     42:   PetscCall(PetscOptionsBool("-view_coord", "View coordinates of element closures", "ex8.c", view_coord, &view_coord, NULL));
     43:   PetscCall(PetscOptionsBoundedInt("-num_fields", "The number of fields to use", "ex8.c", Nf, &Nf, NULL, 1));
     44:   PetscOptionsEnd();

     46:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
     47:   PetscCall(DMSetType(dm, DMPLEX));
     48:   PetscCall(DMSetFromOptions(dm));
     49:   if (project) {
     50:     PetscFE  fe_coords;
     51:     PetscInt cdim;
     52:     PetscCall(DMGetCoordinateDim(dm, &cdim));
     53:     PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, cdim, cdim, PETSC_FALSE, 1, 1, &fe_coords));
     54:     PetscCall(DMSetCoordinateDisc(dm, fe_coords, PETSC_TRUE));
     55:     PetscCall(PetscFEDestroy(&fe_coords));
     56:   }
     57:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
     58:   PetscCall(DMGetDimension(dm, &dim));

     60:   if (Nf == 1) {
     61:     PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, PETSC_DETERMINE, &fe));
     62:     PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
     63:     PetscCall(PetscFEDestroy(&fe));
     64:   } else {
     65:     for (PetscInt f = 0; f < Nf; ++f) {
     66:       char prefix[16];
     67:       PetscCall(PetscSNPrintf(prefix, 16, "f%" PetscInt_FMT "_", f));
     68:       PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, prefix, PETSC_DETERMINE, &fe));
     69:       PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
     70:       PetscCall(PetscFEDestroy(&fe));
     71:     }
     72:   }
     73:   PetscCall(DMCreateDS(dm));
     74:   if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
     75:   PetscCall(DMGetLocalSection(dm, &section));
     76:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
     77:   for (c = cStart; c < cEnd; c++) {
     78:     PetscInt numindices, *indices;
     79:     PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL));
     80:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT "\n", c - cStart));
     81:     PetscCall(PetscIntView(numindices, indices, PETSC_VIEWER_STDOUT_SELF));
     82:     PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL));
     83:   }
     84:   {
     85:     DMLabel         domain_label;
     86:     IS              valueIS;
     87:     const PetscInt *values;
     88:     PetscInt        Nv;

     90:     PetscCall(DMGetLabel(dm, "Face Sets", &domain_label));
     91:     PetscCall(DMLabelGetNumValues(domain_label, &Nv));
     92:     // Check that discretization bas any values on faces
     93:     {
     94:       PetscDS         ds;
     95:       PetscFE         fe;
     96:       IS              fieldIS;
     97:       const PetscInt *fields;
     98:       PetscInt        Nf, dmf = 0, dsf = -1;

    100:       PetscCall(DMGetRegionDS(dm, domain_label, &fieldIS, &ds, NULL));
    101:       // Translate DM field number to DS field number
    102:       PetscCall(ISGetIndices(fieldIS, &fields));
    103:       PetscCall(ISGetSize(fieldIS, &Nf));
    104:       for (PetscInt f = 0; f < Nf; ++f) {
    105:         if (dmf == fields[f]) {
    106:           dsf = f;
    107:           break;
    108:         }
    109:       }
    110:       PetscCall(ISRestoreIndices(fieldIS, &fields));
    111:       PetscCall(PetscDSGetDiscretization(ds, dsf, (PetscObject *)&fe));
    112:       PetscCall(PetscFEGetHeightSubspace(fe, 1, &fe));
    113:       if (!fe) Nv = 0;
    114:     }
    115:     PetscCall(DMLabelGetValueIS(domain_label, &valueIS));
    116:     PetscCall(ISGetIndices(valueIS, &values));
    117:     for (PetscInt v = 0; v < Nv; ++v) {
    118:       PetscInt *elem_restr_offsets_face;
    119:       PetscInt  num_elem, elem_size, num_comp, num_dof;

    121:       PetscCall(DMPlexGetLocalOffsets(dm, domain_label, values[v], 1, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets_face));
    122:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "========= Face restriction; marker: %" PetscInt_FMT " ========\n", values[v]));
    123:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "num_elem: %" PetscInt_FMT ", elem_size: %" PetscInt_FMT ", num_dof:  %" PetscInt_FMT "\n", num_elem, elem_size, num_dof));
    124:       for (PetscInt c = 0; c < num_elem; c++) PetscCall(PetscIntView(elem_size, &elem_restr_offsets_face[c * elem_size], PETSC_VIEWER_STDOUT_SELF));
    125:       PetscCall(PetscFree(elem_restr_offsets_face));
    126:     }
    127:     PetscCall(ISRestoreIndices(valueIS, &values));
    128:     PetscCall(ISDestroy(&valueIS));
    129:   }
    130:   if (view_coord) {
    131:     DM        cdm, cell_dm;
    132:     Vec       X;
    133:     PetscInt  cdim;
    134:     PetscBool sparseLocalize;

    136:     PetscCall(DMGetCoordinatesLocalSetUp(dm));
    137:     PetscCall(DMGetCoordinateDim(dm, &cdim));
    138:     PetscCall(DMGetCoordinateDM(dm, &cdm));
    139:     PetscCall(PetscObjectSetName((PetscObject)cdm, "coords"));
    140:     if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
    141:     PetscCall(DMLocalizeCoordinates(dm));
    142:     PetscCall(DMGetSparseLocalize(dm, &sparseLocalize));
    143:     for (c = cStart; c < cEnd; ++c) {
    144:       const PetscScalar *array;
    145:       PetscScalar       *x = NULL;
    146:       PetscInt           ndof;
    147:       PetscBool          isDG;

    149:       PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
    150:       PetscCheck(ndof % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "ndof not divisible by cdim");
    151:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT " coordinates\n", c - cStart));
    152:       for (PetscInt i = 0; i < ndof; i += cdim) PetscCall(PetscScalarView(cdim, &x[i], PETSC_VIEWER_STDOUT_SELF));
    153:       PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
    154:     }
    155:     PetscCall(ViewOffsets(dm, NULL));
    156:     PetscCall(DMGetCoordinatesLocal(dm, &X));
    157:     PetscCall(ViewOffsets(cdm, X));
    158:     PetscCall(DMGetCellCoordinateDM(dm, &cell_dm));
    159:     PetscCall(PetscObjectSetName((PetscObject)cell_dm, "cell coords"));
    160:     if (cell_dm && !sparseLocalize) {
    161:       PetscCall(DMGetCellCoordinatesLocal(dm, &X));
    162:       PetscCall(ViewOffsets(cell_dm, X));
    163:     }
    164:   }
    165:   PetscCall(DMDestroy(&dm));
    166:   PetscCall(PetscFinalize());
    167:   return 0;
    168: }

    170: /*TEST

    172:   test:
    173:     suffix: 1d_q2
    174:     args: -dm_plex_dim 1 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2
    175:   test:
    176:     suffix: 2d_q1
    177:     args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
    178:   test:
    179:     suffix: 2d_q1d
    180:     args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -petscdualspace_lagrange_continuity 0
    181:   test:
    182:     suffix: 2d_q1_q1d
    183:     args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 1 -f1_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -f1_petscdualspace_lagrange_continuity 0
    184:   test:
    185:     suffix: 2d_q2
    186:     args: -dm_plex_dim 2 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
    187:   test:
    188:     suffix: 2d_q2_q1
    189:     args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 2 -f1_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,1
    190:   test:
    191:     suffix: 2d_q2_p0
    192:     args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 2 -f1_petscspace_degree 0 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -f1_petscdualspace_lagrange_continuity 0
    193:   test:
    194:     suffix: 2d_q2_q1d
    195:     args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -num_fields 2 \
    196:             -f0_petscspace_degree 2 \
    197:             -f1_petscspace_degree 1 -f1_petscdualspace_lagrange_continuity 0
    198:   test:
    199:     suffix: 2d_q2_p1d
    200:     args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -num_fields 2 \
    201:             -f0_petscspace_degree 2 \
    202:             -f1_petscspace_degree 1 -f1_petscspace_poly_tensor 0 -f1_petscdualspace_lagrange_continuity 0
    203:   test:
    204:     suffix: 2d_q3
    205:     args: -dm_plex_dim 2 -petscspace_degree 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1

    207:   test:
    208:     suffix: 3d_q1
    209:     args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
    210:   test:
    211:     suffix: 1d_q1_periodic
    212:     requires: !complex
    213:     args: -dm_plex_dim 1 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3 -dm_plex_box_bd periodic -dm_view -view_coord
    214:   testset:
    215:     suffix: 2d_q1_periodic
    216:     requires: !complex
    217:     args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 3,2 -dm_plex_box_bd periodic,none \
    218:           -petscspace_degree 1 -dm_view -view_coord

    220:     test:
    221:     test:
    222:       suffix: sparse
    223:       args: -dm_sparse_localize false -dm_localize 0

    225:   test:
    226:     suffix: 3d_q1_periodic
    227:     requires: !complex
    228:     args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3,2,1 -dm_plex_box_bd periodic,none,none -dm_view -view_coord
    229:   test:
    230:     suffix: 3d_q1_periodic_project
    231:     requires: !complex
    232:     args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,3 -dm_plex_box_bd none,none,periodic -dm_view -view_coord -project_coordinates

    234:   test:
    235:     suffix: 3d_q2_periodic # not actually periodic because only 2 cells
    236:     args: -dm_plex_dim 3 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2,2 -dm_plex_box_bd periodic,none,periodic -dm_view

    238: TEST*/