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