1 static char help[] = "Test FEM layout with constraints\n\n"; 2 3 #include <petsc.h> 4 5 int main(int argc,char **argv) { 6 DM dm,pdm; 7 PetscSection section; 8 const PetscInt field = 0; 9 char ifilename[PETSC_MAX_PATH_LEN]; 10 PetscInt sdim,s,pStart,pEnd,p,numVS,numPoints; 11 PetscInt constraints[1]; 12 IS setIS,pointIS; 13 const PetscInt *setID,*pointID; 14 15 PetscFunctionBeginUser; 16 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 17 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"FEM Layout Options","ex98"); 18 PetscCall(PetscOptionsString("-i","Filename to read","ex98",ifilename,ifilename,sizeof(ifilename),NULL)); 19 PetscOptionsEnd(); 20 21 PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD,ifilename,NULL,PETSC_TRUE,&dm)); 22 PetscCall(DMPlexDistributeSetDefault(dm,PETSC_FALSE)); 23 PetscCall(DMSetFromOptions(dm)); 24 25 PetscCall(DMPlexDistribute(dm,0,NULL,&pdm)); 26 if (pdm) { 27 PetscCall(DMDestroy(&dm)); 28 dm = pdm; 29 } 30 PetscCall(DMViewFromOptions(dm,NULL,"-dm_view")); 31 32 /* create a section */ 33 PetscCall(DMGetDimension(dm,&sdim)); 34 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject) dm),§ion)); 35 PetscCall(PetscSectionSetNumFields(section,1)); 36 PetscCall(PetscSectionSetFieldName(section,field,"U")); 37 PetscCall(PetscSectionSetFieldComponents(section,field,sdim)); 38 PetscCall(DMPlexGetChart(dm,&pStart,&pEnd)); 39 PetscCall(PetscSectionSetChart(section,pStart,pEnd)); 40 41 /* initialize the section storage for a P1 field */ 42 PetscCall(DMPlexGetDepthStratum(dm,0,&pStart,&pEnd)); 43 for (p = pStart; p < pEnd; ++p) { 44 PetscCall(PetscSectionSetDof(section,p,sdim)); 45 PetscCall(PetscSectionSetFieldDof(section,p,0,sdim)); 46 } 47 48 /* add constraints at all vertices belonging to a vertex set */ 49 /* first pass is to reserve space */ 50 PetscCall(DMGetLabelSize(dm,"Vertex Sets",&numVS)); 51 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"# Vertex set: %" PetscInt_FMT "\n",numVS)); 52 PetscCall(DMGetLabelIdIS(dm,"Vertex Sets",&setIS)); 53 PetscCall(ISGetIndices(setIS,&setID)); 54 for (s = 0; s < numVS; ++s){ 55 PetscCall(DMGetStratumIS(dm,"Vertex Sets",setID[s],&pointIS)); 56 PetscCall(DMGetStratumSize(dm,"Vertex Sets",setID[s],&numPoints)); 57 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"set %" PetscInt_FMT " size: %" PetscInt_FMT "\n",s,numPoints)); 58 PetscCall(ISGetIndices(pointIS,&pointID)); 59 for (p = 0; p < numPoints; ++p){ 60 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\t point %" PetscInt_FMT "\n",pointID[p])); 61 PetscCall(PetscSectionSetConstraintDof(section,pointID[p],1)); 62 PetscCall(PetscSectionSetFieldConstraintDof(section,pointID[p],field,1)); 63 } 64 PetscCall(ISRestoreIndices(pointIS,&pointID)); 65 PetscCall(ISDestroy(&pointIS)); 66 } 67 68 PetscCall(PetscSectionSetUp(section)); 69 70 /* add constraints at all vertices belonging to a vertex set */ 71 /* second pass is to assign constraints to a specific component / dof */ 72 for (s = 0; s < numVS; ++s){ 73 PetscCall(DMGetStratumIS(dm,"Vertex Sets",setID[s],&pointIS)); 74 PetscCall(DMGetStratumSize(dm,"Vertex Sets",setID[s],&numPoints)); 75 PetscCall(ISGetIndices(pointIS,&pointID)); 76 for (p = 0; p < numPoints; ++p){ 77 constraints[0] = setID[s]%sdim; 78 PetscCall(PetscSectionSetConstraintIndices(section,pointID[p],constraints)); 79 PetscCall(PetscSectionSetFieldConstraintIndices(section,pointID[p],field,constraints)); 80 } 81 PetscCall(ISRestoreIndices(pointIS,&pointID)); 82 PetscCall(ISDestroy(&pointIS)); 83 } 84 PetscCall(ISRestoreIndices(setIS,&setID)); 85 PetscCall(ISDestroy(&setIS)); 86 PetscCall(PetscObjectViewFromOptions((PetscObject) section,NULL,"-dm_section_view")); 87 88 PetscCall(PetscSectionDestroy(§ion)); 89 PetscCall(DMDestroy(&dm)); 90 91 PetscCall(PetscFinalize()); 92 return 0; 93 } 94 95 /*TEST 96 build: 97 requires: exodusii pnetcdf !complex 98 testset: 99 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/SquareFaceSet.exo -dm_view -dm_section_view 100 nsize: 1 101 102 test: 103 suffix: 0 104 args: 105 106 TEST*/ 107