1 static char help[] = "Element closure restrictions in tensor/lexicographic/spectral-element ordering using DMPlex\n\n"; 2 3 #include <petscdmplex.h> 4 5 int main(int argc, char **argv) 6 { 7 DM dm; 8 PetscSection section; 9 PetscFE fe; 10 PetscInt dim,c,cStart,cEnd; 11 PetscBool view_coord = PETSC_FALSE; 12 PetscErrorCode ierr; 13 14 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 15 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Tensor closure restrictions", "DMPLEX");CHKERRQ(ierr); 16 ierr = PetscOptionsBool("-view_coord", "View coordinates of element closures", "ex8.c", view_coord, &view_coord, NULL);CHKERRQ(ierr); 17 ierr = PetscOptionsEnd();CHKERRQ(ierr); 18 19 ierr = DMCreate(PETSC_COMM_WORLD, &dm);CHKERRQ(ierr); 20 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 21 ierr = DMSetFromOptions(dm);CHKERRQ(ierr); 22 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr); 23 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 24 25 ierr = PetscFECreateDefault(PETSC_COMM_SELF,dim,1,PETSC_FALSE,NULL,PETSC_DETERMINE,&fe);CHKERRQ(ierr); 26 ierr = DMAddField(dm,NULL,(PetscObject)fe);CHKERRQ(ierr); 27 ierr = DMCreateDS(dm);CHKERRQ(ierr); 28 ierr = DMPlexSetClosurePermutationTensor(dm,PETSC_DETERMINE,NULL);CHKERRQ(ierr); 29 ierr = DMGetLocalSection(dm,§ion);CHKERRQ(ierr); 30 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 31 for (c=cStart; c<cEnd; c++) { 32 PetscInt numindices,*indices; 33 ierr = DMPlexGetClosureIndices(dm,section,section,c,PETSC_TRUE,&numindices,&indices,NULL,NULL);CHKERRQ(ierr); 34 ierr = PetscPrintf(PETSC_COMM_SELF,"Element #%" PetscInt_FMT "\n",c-cStart);CHKERRQ(ierr); 35 ierr = PetscIntView(numindices,indices,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 36 ierr = DMPlexRestoreClosureIndices(dm,section,section,c,PETSC_TRUE,&numindices,&indices,NULL,NULL);CHKERRQ(ierr); 37 } 38 if (view_coord) { 39 DM cdm; 40 Vec X; 41 PetscInt cdim; 42 ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr); 43 ierr = DMPlexSetClosurePermutationTensor(cdm,PETSC_DETERMINE,NULL);CHKERRQ(ierr); 44 ierr = DMGetCoordinatesLocal(dm, &X);CHKERRQ(ierr); 45 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 46 for (c=cStart; c<cEnd; c++) { 47 PetscScalar *x; 48 PetscInt ndof; 49 ierr = DMPlexVecGetClosure(cdm, NULL, X, c, &ndof, &x);CHKERRQ(ierr); 50 PetscCheck(ndof % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "ndof not divisible by cdim"); 51 ierr = PetscPrintf(PETSC_COMM_SELF,"Element #%" PetscInt_FMT " coordinates\n",c-cStart);CHKERRQ(ierr); 52 for (PetscInt i=0; i<ndof; i+= cdim) { 53 ierr = PetscScalarView(cdim, &x[i], PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 54 } 55 ierr = DMPlexVecRestoreClosure(cdm, NULL, X, c, &ndof, &x);CHKERRQ(ierr); 56 } 57 } 58 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 59 ierr = DMDestroy(&dm);CHKERRQ(ierr); 60 ierr = PetscFinalize(); 61 return ierr; 62 } 63 64 /*TEST 65 66 test: 67 suffix: 1d_q2 68 args: -dm_plex_dim 1 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2 69 test: 70 suffix: 2d_q1 71 args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 72 test: 73 suffix: 2d_q2 74 args: -dm_plex_dim 2 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 75 test: 76 suffix: 2d_q3 77 args: -dm_plex_dim 2 -petscspace_degree 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 78 test: 79 suffix: 3d_q1 80 args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 81 test: 82 suffix: 1d_q1_periodic 83 args: -dm_plex_dim 1 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2 -dm_plex_box_bd periodic -dm_view 84 test: 85 suffix: 2d_q1_periodic 86 args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -dm_plex_box_bd periodic,none -dm_view 87 test: 88 suffix: 3d_q2_periodic 89 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 90 91 TEST*/ 92