1 static char help[] = "Element closure restrictions in tensor/lexicographic/spectral-element ordering using DMPlex\n\n"; 2 3 #include <petscdmplex.h> 4 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; 11 12 PetscFunctionBegin; 13 PetscCall(PetscObjectGetName((PetscObject)dm, &name)); 14 PetscCall(DMPlexGetLocalOffsets(dm, NULL, 0, 0, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets)); 15 PetscCall(PetscPrintf(PETSC_COMM_SELF,"DM %s offsets: num_elem %" PetscInt_FMT ", size %" PetscInt_FMT 16 ", 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++) { 22 PetscCall(PetscScalarView(num_comp, &x[elem_restr_offsets[c*elem_size+i]], PETSC_VIEWER_STDERR_SELF)); 23 } 24 } 25 } 26 if (X) PetscCall(VecRestoreArrayRead(X, &x)); 27 PetscCall(PetscFree(elem_restr_offsets)); 28 PetscFunctionReturn(0); 29 } 30 31 int main(int argc, char **argv) 32 { 33 DM dm; 34 PetscSection section; 35 PetscFE fe; 36 PetscInt dim,c,cStart,cEnd; 37 PetscBool view_coord = PETSC_FALSE, tensor = PETSC_TRUE, project = PETSC_FALSE; 38 39 PetscFunctionBeginUser; 40 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 41 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Tensor closure restrictions", "DMPLEX"); 42 PetscCall(PetscOptionsBool("-closure_tensor", "Apply DMPlexSetClosurePermutationTensor", "ex8.c", tensor, &tensor, NULL)); 43 PetscCall(PetscOptionsBool("-project_coordinates", "Call DMProjectCoordinates explicitly", "ex8.c", project, &project, NULL)); 44 PetscCall(PetscOptionsBool("-view_coord", "View coordinates of element closures", "ex8.c", view_coord, &view_coord, NULL)); 45 PetscOptionsEnd(); 46 47 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 48 PetscCall(DMSetType(dm, DMPLEX)); 49 PetscCall(DMSetFromOptions(dm)); 50 if (project) { 51 PetscFE fe_coords; 52 PetscInt cdim; 53 PetscCall(DMGetCoordinateDim(dm, &cdim)); 54 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, cdim, cdim, PETSC_FALSE, 1, 1, &fe_coords)); 55 PetscCall(DMProjectCoordinates(dm, fe_coords)); 56 PetscCall(PetscFEDestroy(&fe_coords)); 57 } 58 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 59 PetscCall(DMGetDimension(dm, &dim)); 60 61 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF,dim,1,PETSC_FALSE,NULL,PETSC_DETERMINE,&fe)); 62 PetscCall(DMAddField(dm,NULL,(PetscObject)fe)); 63 PetscCall(DMCreateDS(dm)); 64 if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm,PETSC_DETERMINE,NULL)); 65 PetscCall(DMGetLocalSection(dm,§ion)); 66 PetscCall(DMPlexGetHeightStratum(dm,0,&cStart,&cEnd)); 67 for (c=cStart; c<cEnd; c++) { 68 PetscInt numindices,*indices; 69 PetscCall(DMPlexGetClosureIndices(dm,section,section,c,PETSC_TRUE,&numindices,&indices,NULL,NULL)); 70 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Element #%" PetscInt_FMT "\n",c-cStart)); 71 PetscCall(PetscIntView(numindices,indices,PETSC_VIEWER_STDOUT_SELF)); 72 PetscCall(DMPlexRestoreClosureIndices(dm,section,section,c,PETSC_TRUE,&numindices,&indices,NULL,NULL)); 73 } 74 if (view_coord) { 75 DM cdm; 76 Vec X; 77 PetscInt cdim; 78 79 PetscCall(DMGetCoordinateDim(dm, &cdim)); 80 PetscCall(DMGetCoordinateDM(dm, &cdm)); 81 PetscCall(PetscObjectSetName((PetscObject) cdm, "coords")); 82 if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL)); 83 for (c = cStart; c < cEnd; ++c) { 84 const PetscScalar *array; 85 PetscScalar *x = NULL; 86 PetscInt ndof; 87 PetscBool isDG; 88 89 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &ndof, &array, &x)); 90 PetscCheck(ndof % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "ndof not divisible by cdim"); 91 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT " coordinates\n", c-cStart)); 92 for (PetscInt i = 0; i < ndof; i+= cdim) PetscCall(PetscScalarView(cdim, &x[i], PETSC_VIEWER_STDOUT_SELF)); 93 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &ndof, &array, &x)); 94 } 95 PetscCall(ViewOffsets(dm, NULL)); 96 PetscCall(DMGetCoordinatesLocal(dm, &X)); 97 PetscCall(ViewOffsets(cdm, X)); 98 } 99 PetscCall(PetscFEDestroy(&fe)); 100 PetscCall(DMDestroy(&dm)); 101 PetscCall(PetscFinalize()); 102 return 0; 103 } 104 105 /*TEST 106 107 test: 108 suffix: 1d_q2 109 args: -dm_plex_dim 1 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2 110 test: 111 suffix: 2d_q1 112 args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 113 test: 114 suffix: 2d_q2 115 args: -dm_plex_dim 2 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 116 test: 117 suffix: 2d_q3 118 args: -dm_plex_dim 2 -petscspace_degree 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 119 test: 120 suffix: 3d_q1 121 args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 122 test: 123 suffix: 1d_q1_periodic 124 requires: !complex 125 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 126 test: 127 suffix: 2d_q1_periodic 128 requires: !complex 129 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 130 test: 131 suffix: 3d_q1_periodic 132 requires: !complex 133 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 134 test: 135 suffix: 3d_q1_periodic_project 136 requires: !complex 137 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 138 139 test: 140 suffix: 3d_q2_periodic # not actually periodic because only 2 cells 141 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 142 143 TEST*/ 144