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

  3: #include <petscdmplex.h>

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

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

 32: int main(int argc, char **argv)
 33: {
 34:   DM             dm;
 35:   PetscSection   section;
 36:   PetscFE        fe;
 37:   PetscInt       dim,c,cStart,cEnd;
 38:   PetscBool      view_coord = PETSC_FALSE, tensor = PETSC_TRUE;

 41:   PetscInitialize(&argc,&argv,NULL,help);
 42:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Tensor closure restrictions", "DMPLEX");
 43:   PetscOptionsBool("-closure_tensor", "Apply DMPlexSetClosurePermutationTensor", "ex8.c", tensor, &tensor, NULL);
 44:   PetscOptionsBool("-view_coord", "View coordinates of element closures", "ex8.c", view_coord, &view_coord, NULL);
 45:   PetscOptionsEnd();

 47:   DMCreate(PETSC_COMM_WORLD, &dm);
 48:   DMSetType(dm, DMPLEX);
 49:   DMSetFromOptions(dm);
 50:   DMViewFromOptions(dm, NULL, "-dm_view");
 51:   DMGetDimension(dm, &dim);

 53:   PetscFECreateDefault(PETSC_COMM_SELF,dim,1,PETSC_FALSE,NULL,PETSC_DETERMINE,&fe);
 54:   DMAddField(dm,NULL,(PetscObject)fe);
 55:   DMCreateDS(dm);
 56:   if (tensor) DMPlexSetClosurePermutationTensor(dm,PETSC_DETERMINE,NULL);
 57:   DMGetLocalSection(dm,&section);
 58:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
 59:   for (c=cStart; c<cEnd; c++) {
 60:     PetscInt numindices,*indices;
 61:     DMPlexGetClosureIndices(dm,section,section,c,PETSC_TRUE,&numindices,&indices,NULL,NULL);
 62:     PetscPrintf(PETSC_COMM_SELF,"Element #%" PetscInt_FMT "\n",c-cStart);
 63:     PetscIntView(numindices,indices,PETSC_VIEWER_STDOUT_SELF);
 64:     DMPlexRestoreClosureIndices(dm,section,section,c,PETSC_TRUE,&numindices,&indices,NULL,NULL);
 65:   }
 66:   if (view_coord) {
 67:     DM cdm;
 68:     Vec X;
 69:     PetscInt cdim;
 70:     DMGetCoordinateDM(dm,&cdm);
 71:     PetscObjectSetName((PetscObject)cdm, "coords");
 72:     if (tensor) DMPlexSetClosurePermutationTensor(cdm,PETSC_DETERMINE,NULL);
 73:     DMGetCoordinatesLocal(dm, &X);
 74:     DMGetCoordinateDim(dm, &cdim);
 75:     for (c=cStart; c<cEnd; c++) {
 76:       PetscScalar *x = NULL;
 77:       PetscInt ndof;
 78:       DMPlexVecGetClosure(cdm, NULL, X, c, &ndof, &x);
 80:       PetscPrintf(PETSC_COMM_SELF,"Element #%" PetscInt_FMT " coordinates\n",c-cStart);
 81:       for (PetscInt i=0; i<ndof; i+= cdim) {
 82:         PetscScalarView(cdim, &x[i], PETSC_VIEWER_STDOUT_SELF);
 83:       }
 84:       DMPlexVecRestoreClosure(cdm, NULL, X, c, &ndof, &x);
 85:     }
 86:     ViewOffsets(dm, NULL);
 87:     ViewOffsets(cdm, X);
 88:   }
 89:   PetscFEDestroy(&fe);
 90:   DMDestroy(&dm);
 91:   PetscFinalize();
 92:   return 0;
 93: }

 95: /*TEST

 97:   test:
 98:     suffix: 1d_q2
 99:     args: -dm_plex_dim 1 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2
100:   test:
101:     suffix: 2d_q1
102:     args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
103:   test:
104:     suffix: 2d_q2
105:     args: -dm_plex_dim 2 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
106:   test:
107:     suffix: 2d_q3
108:     args: -dm_plex_dim 2 -petscspace_degree 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1
109:   test:
110:     suffix: 3d_q1
111:     args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
112:   test:
113:     suffix: 1d_q1_periodic
114:     requires: !complex
115:     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
116:   test:
117:     suffix: 2d_q1_periodic
118:     requires: !complex
119:     args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3,2 -dm_plex_box_bd periodic,none -dm_view -view_coord
120:   test:
121:     suffix: 3d_q1_periodic
122:     requires: !complex
123:     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
124:   test:
125:     suffix: 3d_q2_periodic  # not actually periodic because only 2 cells
126:     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

128: TEST*/