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